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SUMMARY 


In  this  study,  \j/e  searched  for  evidence  of  a  strange  attractor  associated 
with  the  saturation  of  middle  atmosphere  gravity  waves  in  the  echo  data  from 
a  partial  reflection  radar  located  in  Saskatoon,  Canada.  Theiler's  extension  of 
the  Grassbcrger-Procaccia  correlation  integral  algorithm  was  used  to  estimate 
the  fractal  dimension  of  the  attractor.  — 

\ 

Chaotic  regimes  have  been  observed  in\  experimental  fluid  studies  of  the 
transition  from  ordered  to  turbulent  behavionWreaking  gravity  waves  are 


thought  to  decay  to  turbulence,  transporting  momentum  from  the  lower  to 
f)1/ 

upper  atmosphere. ^/Extending  Ipe— re$ults~frem?>  laboratory  h  Studies^*  to  middle 
atmosphere  gravity  waves,  it  seems  reasonable  to  expect  to  find  a  strange 
attractor  in  gravity  wave  saturation. 

Echo  data  was  analyzed  because  it  offered  a  high  sampling  rate.  The 
Grassbcrger-Procaccia  algorithm  places  stringent  requirements  upon  the 
amount  of  data  necessary  to  obtain  an  accurate  estimate  of  the  system 
dimension;  a  large  number  of  points  is  required. 

We  did  not  detect  a  strange  attractor  with  dimension  <3  in  the  data  from 
the  Saskatoon  partial  reflection  radar  for  the  time  scales  (6  min  39  s)  which 

were  studied.-*  Because  of  the  small  number  of  points  which  were  examined.we 

I 

can  not  assert  that  there  was  only  noise  in  the  data.  However,  the  supporting 
evidence  from  the  power  spectra  suggest  that  we  mainly  investigated  time 
scales^,in  the  viscous  and  inertial  regions. 

6-  This  study  can  not  assert  that  a  strange  attractors  is  absent  in  gravity 
wave  absorption.  ^Fhe'^data  requirements  to  implement  the  Grassberger-" 
Procaccia  algorithm  make  it  unlikely  that  such  an  attractor,  if  it  exists,  will  be 
detected.  (Calculations  of  the  amount  of  data  necessary  to  estimate  the 
dimension  ‘indicate  that  over  6  hours  of  data  would  be  required  to  detect  a 

t 

I 

strange  attractor  in  gravity  wave  absorption. 
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CHAPTER  I 


INTRODUCTION 


In  the  last  10  years,  MLT  (mesosphere,  lower  thermosphere)  radars 
have  been  increasingly  important  in  the  observation  of  the  middle 
atmosphere1.  During  this  same  period,  gravity  waves  and  gravity  wave 
saturation  have  been  recognized  as  playing  a  vital  role  in  the  maintenance 
and  modification  of  the  mesospheric  circulation  and  temperature  distribution 
(Fritts,  1984).  In  addition,  this  same  period  also  saw  great  advances  in  the  area 
popularly  referred  to  as  "chaos  theory". 

Chaotic  behavior  (or  the  presence  of  strange  attractors2)  has  been 
observed  in  experimental  studies  of  the  transition  from  laminar  to  turbulent 
fluid  flow  (see  Swinney  (1983)  for  a  review  of  early  experimental  results  for 
different  systems).  Gravity  waves  undergoing  saturation3  break  down;  the  end 
product  is  turbulence.  Techniques  have  been  developed  to  determine  the 
presence  and  dimension  of  a  strange  attractor  in  a  set  of  data  (Grassbcrger  and 
Procaccia,  1983,  1984).  It  does  not  seem  unreasonable  to  expect  to  find  a  strange 
attractor  associated  with  the  saturation  of  middle  atmosphere  gravity  waves. 

Saturatio  .  of  middle  atmosphere  gravity  waves  has  been  inferred  from 
theory  and  :  ^direct  observations.  Gravity  wave  saturation  currently  provides 
the  only  known  mechanism  for  the  observed  structure  of  the  mesospheric 

1  In  this  thesis,  the  commonly  accepted  definition  of  the  middle  atmosphere, 
i.e.,  that  region  of  the  atmosphere  which  encompasses  the  stratosphere  and 
mesosphere,  will  be  used.  This  definition  roughly  includes  the  region  of  the 
atmosphere-  from  10  to  100  km. 

2  An  attractor  is  defined  as  "strange"  if  its  phase  space  trajectories  diverge 
exponentially  on  the  average. 

3  The  term  "gravity  wave  saturation"  refers  to  any  process  that  acts  to  limit  or 
maintain  constant  wave  amplitudes  with  height. 
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circulation  (most  :mportantly,  the  closing  of  the  mesospheric  jets)  and  the 
observed  temperature  structure.  Examples  of  the  indirect  evidence  of  gravity 
wave  saturation  include  observations  of  regions  of  enhanced  turbulent 
diffusion,  measurements  of  super-adiabatic  lapse  rates  in  the  mesosphere  and 
lower  thermosphere  and  measurements  of  momentum  drag  in  the  mesosphere 
(see  Fritts  (1984)  for  a  detailed  description  and  list  of  references  for  these 
different  observations).  The  most  direct  evidence  comes  from  the  observations 
of  Kelvin-Helmholtz  billows  near  the  summer  mesopause,  as  revealed  in  the 
perturbation  of  noctilucent  clouds. 

Chaos  theory  has  provided  some  insights  into  the  transition  of  fluid 
flows  from  orderly  to  turbulent  regimes.  Traditional  analysis  of  meteorological 
data  has  centered  on  the  search  for  wavelike  or  periodic  behavior.  Tools  such 
as  Fourier  analysis  yield  no  meaningful  results  when  applied  to  aperiodic 
signals.  Irregular  or  aperiodic  signals  usually  are  filtered  out  or  deemed  noise. 
This  noise  may  hide  the  presence  of  a  strange  attractor  within  the  data.  As 
Froehling  ct  al,  (1981)  point  out,  "power  spectral  analysis,  for  example, 
characterizes  aperiodic  behavior  by  the  presence  of  broadband  noise  in  the 
power  spectrum,  but  broadband  noise  can  be  produced  by  systems  requiring 
either  a  small  or  large  number  of  phase  space  dimensions."4 

The  technique  of  Grassberger  and  Procaccia  (1983,  1984)  has  been  used 
to  examine  various  experimental  data  sets.  Atmanspacher  et  al.  (1988) 
employed  the  correlation  integral  technique  of  Grassberger  and  Procaccia  to 
examine  the  chaotic  attractor  associated  with  X-ray  counts  from  the  neutron 
star  Her  X-l.  Tsonis  and  Eisner  (1988)  employed  this  same  technique  on  daytime 
vertical  wind  velocities  in  the  boundary  layer.  BrandstSter  and  Swinney 
(1987)  applied  the  Grassberger-Procaccia  algorithm  to  experimental  data 
obtained  from  the  observation  of  Couette-Taylor  flow.  Elgar  and  Mayer-Kress 
(1989)  applied  the  correlation  integral  technique  to  ocean  wave  data  and  found 


4  Froehling  et  al.  (1981),  p605. 
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a  system  with  a  correlation  dimension^  greater  than  9.  Osborne  et  al.  (1986) 
used  the  Grassberger-Procaccia  algorithm  to  find  a  correlation  dimension  of 
1.4  in  the  motions  of  buoys  in  the  Pacific  ocean.  There  are  many  additional 
studies  in  which  this  technique  was  used;  some  are  described  in  later  chapters. 

l.i  Motivation 

This  study  proposes  to  search  for  a  strange  attractor  associated  with  the 
saturation  of  middle  atmosphere  gravity  waves  in  the  echo  data  obtained  from 
a  partial  reflection  radar  located  in  Saskatoon,  Canada.  This  study  rests  on  the 
hypothesis  that  there  is  a  strange  attractor  associated  with  the  saturation  of 
middle  atmosphere  gravity  waves  and  that  it  can  be  detected  in  middle 
atmosphere  data. 

We  chose  to  use  the  raw  echo  data  from  the  partial  reflection  radar 
because  it  had  a  high  sampling  rate;  the  normal  post-processing  which 
retrieves  horizontal  winds  yields  only  one  data  point  for  approximately  two 
minutes  of  echo  data  and  introduces  a  degree  of  smoothing  to  the  signal.  The 
relationship  between  the  echoes  and  the  dynamics  and  physics  of  the  middle 
atmosphere  is  not  completely  understood;  this  drawback  will  hinder  the 
interpretation  of  the  physical  meaning  of  any  attractor  that  might  be  found. 
The  data  set  studied  here  was  chosen  because  of  the  possibility  that  it  contained 
gravity  waves  and  it  was  given  to  us  (free!). 

The  laboratory  studies  described  in  Chapter  2  suggest  that  there  is  a 
transitional  regime  in  many  fluids  between  laminar  and  random  behavior.. 

This  transitional  regime  occurs  when  some  critical  stability  parameter  is 
exceeded.  Or.ce  this  critical  threshold  is  passed,  the  flow  is  considered  chaotic 
and  is  characterized  by  a  low  dimension,  non-integer  attractor.  The  dimension 
of  these  systems  is  integer  for  stability  parameter  values  which  are  below  the 
critical  level  but  becomes  non-integer  once  this  critical  threshold  is  passed 

5  See  Chapter  III  for  a  discussion  of  'he  different  definitions  of  dimension. 
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and  increases  as  the  stability  of  the  system  decreases.  However,  these  systems 
take  many  different  routes  from  ordered  (i.e.,  periodic)  to  chaotic  behavior; 
there  does  not  appear  to  be  a  universal  route  to  chaotic  behavior.  This  will  be 
discussed  in  detail  in  Chapter  III. 

Upward  propagating  gravity  waves  are  thought  to  grow  until  they 
become  unstable  and  decay  to  turbulence,  transporting  momentum  from  the 
lower  atmosphere  to  the  upper  atmosphere.  Their  stability  changes 
continuously  as  the  gravity  waves  propagate  upward  and  is  a  function  of  the 
wave,  the  atmosphere  through  which  it  propagates  and  the  interaction 
between  the  wave  (or  waves,  as  is  the  more  likely  scenario)  and  the 
atmosphere.  If  the  results  from  laboratory  studies  can  be  extended  to  gravity 
waves  in  the  atmosphere,  the  saturation  of  middle  atmosphere  gravity  waves 
may  also  be  characterized  by  a  transitional  regime  and  hence  a  strange 
attractor.  The  dimension  of  this  attractor  should  be  a  function  of  altitude,  since 
the  stability  of  the  gravity  wave  is  a  function  of  atmospheric  variables  which 
vary  with  height.  The  dimension  of  the  gravity  wave  should  become  non¬ 
integer  once  it  becomes  saturated  and  should  increase  as  the  wave  propagates 
upward  past  the  saturation  level. 

This  study  does  not  depend  on  the  exact  mechanism  of  gravity  wave 
saturation;  there  are  many  different  conceptual  models  of  gravity  wave 
breaking.  It  assumes  only  that  gravity  waves  do  become  saturated  in  the 
middle  atmosphere.  There  are  some  mechanisms  which  limit  gravity  wave 
growth  (e.g.,  nonlinear  wave-wave  interaction)  but  do  not  result  in  wave 
breaking.  In  addition,  other  types  of  fluid  instabilities  occur  in  the  middle 
atmosphere  besides  those  associated  with  gravity  waves,  e.g.,  Kelvin-Helmholtz 
instability.  The  Kelvin-Helmholtz  instability  is  supported  by  observations  of 
billow  clouds  near  the  summer  mesopause  (Fritts,  1984). 

We  focus  on  gravity  waves  because  the  vertical  profiles  generated  by 
the  radar  can  be  used  to  monitor  the  changes  in  the  system  dimension  as  the 
wave  propagates  upward  and  the  stability  changes.  This  does  not  rule  out 
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detecting  strange  attractors  in  the  data  which  are  the  result  of  some  other  type 
of  wave  and  instability.  If  chaotic  behavior  in  the  transition  from  laminar  to 
turbulent  flow  is  universal  then  we  should  find  evidence  of  strange  attractors 
for  fluid  instabilities  other  than  gravity  waves.  Two  different  models  of 
gravity  wave  instability  are  reviewed  in  the  hope  that  finding  a  strange 
attractor  may  offer  insight  into  which  model  better  describes  gravity  wave 
saturation. 

The  hypothesis  that  there  is  a  strange  attractor  associated  with  the 
saturation  of  middle  atmosphere  gravity  waves  which  we  will  be  able  to  detect 
rests  on  many  assumptions.  The  biggest  assumption  is  that  the  transition  to 
chaotic  behavior  observed  in  experimental  studies  is  applicable  to  the  types  of 
fluid  instabilities  that  occur  in  the  atmosphere.  The  second  assumption  is  that 
we  will  be  able  to  detect  the  presence  of  a  strange  attractor  in  atmospheric 
data. 

None  of  the  fluid  studies  that  have  shown  the  transition  from  order  to 
chaos  in  the  laboratory  are  a  particularly  apt  analogy  for  atmospheric  gravity 
waves.  Closed  systems  (e.g.,  Rayleigh-BSnard  convection,  Couette-Taylor  flow) 
are  very  dissimilar  to  gravity  waves.  The  experimental  system  closest  to 
gravity  waves  in  which  the  transition  to  chaotic  behavior  is  observed  is  the 
excited  jet  (Bonetti  and  Boon,  1989).  Chaotic  behavior  may  be  specific  to  these 
systems  and  not  indicative  of  a  more  universal  behavior. 

We  might  not  be  able  to  detect  a  strange  attractor  in  gravity  wave 
saturation  even  if  it  exists.  Experimental  studies  offer  (tic  opportunity  to  make 
a  long  series  of  observations  of  a  fluid  under  precisely  controlled  conditions. 
Such  controlled  conditions  do  not  exist  in  the  atmosphere;  the  atmosphere 
changes  continuously.  The  search  for  strange  attractors  in  atmospheric  data 
has  been  largely  unsuccessful  despite  studies  which  claim  to  find  them  (see 
Chapter  Ii  for  a  summary  of  these  studies  and  our  critique).  Measurements  of 
the  atmosphere  are  rarely  stationary  (in  the  statistical  sense)  and  never 
contain  as  much  data  as  one  would  like.  These  two  problems  create  an  almost 
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insurmountable  obstacle  to  detecting  strange  attractors  in  atmospheric  data  as 
will  be  discussed  in  later  chapters. 

The  data  set  may  not  contain  any  evidence  of  gravity  wave  saturation. 

In  laboratory  studies,  measurements  can  be  made  at  several  locations  in  the 
fluid  (especially  for  closed  systems)  and  over  a  long  period  of  time  under 
precisely  controlled  conditions.  In  the  case  of  the  excited  jet,  measurements 
were  made  at  one  location  as  the  flow  streamed  by.  We  can  not  measure  gravity 
waves  in  a  similar  manner.  Under  ideal  conditions,  gravity  waves  will 
propagate  through  the  volume  of  the  atmosphere  that  is  probed  by  the  radar 
and  thereby  be  detected.  It  may  be  unlikely  that  gravity  waves  will  undergo 
saturation  in  the  volume  of  the  atmosphere  that  is  directly  being  measured. 
However,  the  turbulence  left  behind  by  gravity  wave  saturation  may  be 
detected  as  it  is  advected  over  the  radar  site  by  the  mean  wind. 

1.2  Organization 


This  work  is  divided  into  six  chapters.  The  first  chapter,  "Introduction", 
is  almost  complete  by  this  point.  The  second  chapter  is  titled  "Gravity  Waves 
and  Chaos  in  Fluids"  followed  by  chapter  III,  "The  Methods  of  Analysis".  The 
data  and  its  source  are  reviewed  in  Chapter  IV,  "Overview  of  the  Data".  The 
fifth  chapter,  "Analysis  and  Interpretation",  contains  the  analysis  of  the  data 
and  describes  its  meaning.  The  final  chapter,  "Conclusion  and 
Recommendations  for  Future  Work",  provides  a  summary  of  the  conclusions 
and  suggestions  for  further  work 

Chapter  II  gives  a  brief  review  of  some  of  experimental  work  done  on 
different  types  of  fluid  flows  in  which  chaotic  behavior  was  observed.  It 
describes  some  of  the  most  common  routes  to  chaotic  behavior  observed  in 
fluid  experiments  and  how  they  are  interpreted.  Chapter  II  also  contains  a 
brief  account  of  attempts  at  detecting  strange  attractors  in  atmospheric  data 
and  the  different  flaws  in  many  of  these  studies.  The  chapter  concludes  with  a 
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very  brief  review  of  gravity  wave  theory;  it  focuses  on  two  of  the  many 
possible  mechanisms  behind  gravity  wave  saturation  in  the  upper 
atmosphere. 

Chapter  III  describes  the  methods  used  to  analyze  the  data.  Because 
analyzing  meteorological  data  for  chaotic  behavior  is  relatively  new,  most  of 
the  chapter  is  devoted  to  a  detailed  overview  of  the  Grassberger-Procaccia 
correlation  integral  algorithm.  The  strengths,  weaknesses  and  limitations  of 
this  algorithm  are  thoroughly  reviewed.  The  chapter  concludes  with  a  brief 
overview  of  the  more  conventional  autocorrelation  and  power  spectrum 
analysis  that  will  be  used  to  supplement  the  analysis  of  system  dimension. 

The  data  used  in  this  analysis  is  described  in  Chapter  IV.  The  theory 
behind  partial  reflection  radar  measurement  techniques  is  briefly  reviewed 
because  it  adds  some  insight  into  how  to  interpret  the  data.  The  data  set 
contained  a  number  of  deficiencies  which  limited  the  scope  of  the  analysis; 
these  are  also  described  in  Chapter  IV. 

Chapter  V  covers  the  implementation  of  the  analysis,  the  results  and 
their  interpretation.  The  first  section  describes  the  implementation  of  the 
Grassberger-Procaccia  algorithm  along  with  the  necessary  supporting 
analysis.  A  sample  of  the  results  of  this  analysis  is  given  in  the  following 
section.  The  chapter  concludes  with  the  interpretation  and  discussion  of  the 
results  of  the  analysis  of  the  system  dimension. 

The  conclusions  and  recommendations  for  future  work  are  given  in 
Chapter  VI.  This  thesis  leaves  many  unanswered  questions  which  provide 
ample  room  for  further  research.  While  no  evidence  of  chaotic  behavior  was 
found  in  this  data  set,  there  still  remains  more  work  to  be  done  on  both 
refining  the  analysis  technique  as  applied  to  atmospheric  data  and  searching 
for  chaotic  behavior  in  the  generation  of  atmospheric  turbulence. 

Appendix  A  contains  the  complete  set  of  graphs  depicting  the  results  of 
the  correlation  integral  algorithm  analysis.  These  are  included  in  their 
entirety  to  fully  document  the  negative  results  of  this  study.  The 


corresponding  slopes  of  the  figures  shown  in  Appendix  A  are  included  in 
Appendix  B.  The  slope  of  the  correlation  integral  should  be  equal  to  the  fractal 
dimension  of  the  attractor  if  there  are  sufficient  points  to  fully  saturate  the 
attractor  in  phase  space. 
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CHAPTER  II 


'  ( 


GRAVITY  WAVES  AND  CHAOS  IN  FLUIDS 

This  chapter  reviews  some  of  the  work  on  chaotic  behavior  in  different 
types  of  fluid  regimes  observed  in  laboratory  experiments.  It  also  provides  a 
brief  description  of  attempts  to  i ind  evidence  of  a  fractal  dimension  (and  thus 
a  strange  attractor)  in  atmospheric  data.  It  concludes  with  a  general 
description  of  internal  gravity  waves  in  the  atmosphere  and  a  discussion  of 
two  possible  mechanisms  of  gravity  wave  breaking. 


Laboratory  experiments  which  examine  chaotic  behavior  in  fluids  can 
be  divided  into  two  categories:  those  in  open  systems  and  those  in  closed 
systems.  In  closed  systems,  the  fluid  is  confined  between  rigid  boundaries.  In 
open  systems,  the  fluid  is  either  not  bound  by  rigid  boundaries  or  the 
boundaries  are  far  enough  away  as  to  not  influence  the  flow.  Couette-Taylor 
flow  and  Rayleigh-Bgnard  convection  arc  examples  of  closed  systems.  The 
excited  jet  is  an  open  system  whose  description  will  follow  that  of  the  closed 
systems. 


Two  of  the  most  frequently  examined  closed  systems  are  Rayleigh- 
Benard  convection  and  Couette-Taylor  How.  Both  systems  provide  well  defined 
examples  of  the  transition  :o  chaotic  behavior  as  exemplified  by  weakly 
developed  turbulence. 


9 


In  the  Rayleigh-Bdnard  system,  the  fluid  is  confined  between  two 
parallel  plates  which  are  held  at  different  temperatures,  usually  by  heating 
the  lower  plate.  The  fluid  develops  convective  cells  whose  behavior  is  a 
function  of  the  dimensionless  Rayleigh  number,  Rd  =  (gad3/Kv)AT,  where  g  is 
the  acceleration  due  to  gravity,  a  the  thermal  expansion  coefficient,  d  the 
distance  between  the  two  plates,  k  the  thermal  diffusivity,  v  the  kinematic 
viscosity  and  AT  the  temperature  boundary  conditions  on  the  side  walls 
(Swinney,  1983). 

In  the  Couette-Taylor  system,  the  fluid  is  confined  between  two 
concentric  cylinders  which  rotate  independently  at  angular  velocity  £2j  and 
£2o-  Most  studies  have  focused  on  the  case  where  the  rotation  rate  of  the  inner 
cylinder  is  zero.  The  behavior  of  Couette-Taylor  flow  is  governed  by  the 
dimensionless  Reynold's  number,  R=((b  -  a)bQo/v).  where  a  and  b  are  the 
radii  of  the  inner  and  outer  cylinders  respectively,  Qo  is  the  angular  velocity 
of  the  outer  cylinder,  and  v  is  the  kinematic  viscosity  (Swinney,  1983). 

2.1. l.l  RaylfcighT-Benard.  Convection 

Bonctti  and  Boon  (1989)  note  that  the  low  dimension  chaotic  attractor  in 
Rayleigh-Benard  convection  is  "a  consequence  of  the  high  confinement 
imposed  by  the  boundaries  an  the  internal  flow  which  results  in  strong 
coupling  between  modes"1.  The  end  result  is  spatially  coherent  "frozen"  flow 
which  is  described  by  a  single,  low  dimension,  chaotic  attractor.  In  ..mil 
aspect  ratio2  Rayleigh-Benard  systems,  the  primary  routes  to  chaos  are  period 
doubling  and  intermittency  (Behringer,  1985);  see  the  following  section  for 
definitions  of  the  different  routes  to  chaos. 

In  large  aspect  ratio  Raylcigh-B6nard  systems,  chaos  is  generated  by 
competition  between  different  unstable  modes,  each  of  which  can  be  described 

1  Bonetti  and  Boon  (1989),  p3322. 

2  The  aspect  ratio  is  defined  as  the  ratio  of  the  horizontal  dimension  to  the 
depth. 
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by  a  localized  chaotic  attractor  (Walden  et  al.,  1985).  The  turbulence  in  these 
systems  is  "spatio-temporal"  which  results  in  a  loss  of  spatial  coherence  in  the 
flow.  Behringer  (1985)  feels  the  precise  route  to  chaos  is  still  unclear. 

Libchaber  et  al.  (1983)  studied  the  route  to  chaos  for  Rayleigh-Bdnard 
convection  in  the  presence  of  a  magnetic  field  as  a  function  of  two  control 
parameters,  the  Rayleigh  number  Rd  and  the  Chandrasekhar  number  Q;  Q  is 
defined  as 


(2.1) 

where  a  is  the  electrical  conductivity,  Bo  is  the  magnitude  of  the  horizontal 
magnetic  field,  d  is  the  depth  of  the  fluid,  p  is  the  fluid  density  and  v  is  the 
kinematic  viscosity.  The  magnetic  field  tends  to  "stiffen"  the  fluid  allowing 
Rayleigh  numbers  higher  than  the  normal  critical  values  to  be  investigated; 
thus,  larger  nonlincaritics  in  the  convection  can  be  examined.  Libchaber  et 
al.  (1983)  found  period  doubling  and  frequency  locking  to  be  the  routes  to 

chaotic  behavior  for  low  Rayleigh  numbers  and  low  magnetic  fields. 

Libchaber  et  al.  (1983)  found  quasi-periodicity  and  soft  mode  instability  (i.e., 

the  interaction  between  oscillatory  instability  and  stationary  instability)  as 
the  routes  to  chaos  for  high  Rayleigh  numbers. 

2.1. 1.2  Couette-Tavlor  Flow 

Brandstater  and  Swinncy  (1987)  examined  chaotic  behavior  in  Couette- 
Taylor  flow.  For  values  of  R/Rc  <  11.7,  where  R  and  Rc  are  the  Reynolds  and 

critical  Reynolds  number  for  the  system,  the  dimension  of  the  system  is  2  for 

modulated,  wavy  vortex  flow.  When  the  ratio  of  Reynolds  numbers  exceeds  that 
threshold,  the  dimension  of  the  system  becomes  non-integer  and  slightly 
greater  than  2.  This  threshold  also  marks  the  first  appearance  of  broadband 
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noise  in  the  power  spectrum.  BrandstSter  and  Swinney  (1987)  note  that  the 
exponential  decay  in  the  power  spectrum  provided  further  evidence  of  non¬ 
periodic  behavior  corresponding  to  a  low  dimension  chaotic  attractor  rather 
than  stochastic  processes.  The  dimension  of  the  system  increased  as  the  ratio 
of  Reynolds  numbers  (R/Rc)  increased  above  the  critical  threshold 
(BrandstSter  et  al.,  1983;  BrandstSter  and  Swinney,  1987). 

BrandstSter  and  Swinney  (1987)  found  that  the  attractor  dimension 
characterized  flow  over  the  entire  annulus.  The  dimension  was  the  same 
(within  the  error  limits  of  the  calculation)  for  measurements  made  at  a 
number  of  different  locations  in  the  flow.  Thus,  the  attractor  characterized  the 
entire  flow  in  the  annulus  rather  a  specific  location  in  'he  fluid.  This  result 
will  not  be  true  of  open  systems. 

Brandstater  and  Swinney  (1987)  found  that  none  of  the  well  established 
routes  to  chaos,  e.g.,  period  doubling,  intermittency,  described  the  transition  to 
chaos  in  this  system  for  the  conditions  they  investigated.  They  speculate  that 
their  experiment  revealed  another  route  to  chaos  which  requires  further 
study. 

2.1.2  Open  Systems 

The  transition  to  fully  developed  turbulence  is  generally  investigated  in 
open  systems  whereas  the  transition  to  weakly  developed  turbulence  is  studied 
in  closed  systems.  The  excited  jet  provides  an  example  of  chaotic  behavior  in 
an  open  system  (Bonetti  and  Boon,  1989).  The  excited  jet  is  of  interest  to  us 
because  of  similarities  to  shear  flows  in  the  atmosphere. 

Bonetti  and  Boon  (1989)  observe  that  the  region  of  growth  of  the  most 
unstable  mode  in  open  systems  is  followed  by  nonlinear  saturation  of  those 
modes  which  generate  advected  coherent  structures.  This  process  leads  to 
three  dimensional  destabilization  and  breakdown  of  these  coherent  structures. 
Bonetti  and  Boon  (1989)  investigated  this  highly  transitional  region  in  open 
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flow  as  exemplified  by  the  excited  jet  in  an  effort  to  examine  the  spatial 
destabilization  of  "coherent  structures"  in  the  flow. 

The  excited  jet  is  a  stream  of  air  that  is  emitted  under  pressure  from  a 
pipe.  The  flow  is  axially  symmetric  and  has  a  Poiseuille  velocity  profile  as  a 
function  of  radial  distance  from  the  center  of  the  flow.  The  jet  is  excited  by 
applying  a  perturbation  to  the  flow  (i.e.,  by  vibrating  the  end  of  the  pipe  from 
which  the  flow  emerges). 

The  excited  jet  has  three  distinct  regions  downstream:  the  laminar 
region,  the  weakly  turbulent  zone  and  the  turbulent  zone.  (Bonetti  and  Boon, 
1989).  The  laminar  region,  nearest  the  source  of  the  flow,  is  characterized  by 
stationary  macroscopic  structures.  This  region  is  followed  by  a  weakly 
turbulent  zone  where  the  macroscopic  structures  are  no  longer  steady  in  time. 
Farther  downstream  is  the  turbulent  zone  where  the  macroscopic  structures 
have  disappeared  and  flow  is  essentially  random:  mixing  occurs  in  this  zone. 
Bonetti  and  Boon  (1989)  note  that  the  appearance  of  these  regions  was 
independent  of  the  excitation  frequency,  although  varying  the  frequency  did 
alter  their  length.  They  sampled  over  8000  periods  of  the  attractor  making 
between  10  and  30  measurements  per  period. 

Although  they  could  not  accurately  determine  the  Kolmogorov  entropy, 
Bonetti  and  Boon  (1989)  did  determine  that  it  had  a  finite  non-zero  value  which 
is  indicative  of  chaotic  behavior  (Grassberger  and  Procaccia,  1983).  Bonetti 
and  Boon  (1989)  found  that  the  flow  was  characterized  by  a  non-integer 
dimension  which  had  an  initial  value  of  less  than  3  but  increased  farther 
downstream  to  between  3  and  4.  The  increase  in  attractor  dimension 
downstream  was  associated  with  a  corresponding  growth  in  broadband  noise 
in  the  power  spectrum.  The  turbulent  region  was  characterized  by  a 
continuous  growth  in  the  correlation  dimension  downstream;  the  attractor  in 
this  region  was  not  saturated  due  to  an  inadequate  number  of  points.  This  was 
most  likely  because  the  flow  became  essentially  random  although  the  limited 
number  of  points  in  the  data  set  makes  this  conclusion  tenuous. 
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Unlike  closed  systems  where  the  attractor  characterized  the  flow  in  the 
entire  system,  the  attractor  dimension  in  open  systems  only  characterized  the 
flow  over  a  local  length  scale  (Bonetti  and  Boon,  1989).  The  local  length  scale  is 
defined  as  the  distance  to  the  first  zero  in  the  spatial  autocorrelation  or  the 
distance  to  the  first  local  minimum  in  the  spatial  mutual  information. 
Intuitively,  this  makes  sense.  The  attractor  dimension  in  closed  systems 
increases  as  an  external  stability  criteria  decreases;  the  change  in  stability 
characterizes  the  entire  system.  The  instability  in  the  open  jet  amplifies 
downstream  corresponding  to  a  continuum  of  stability  changes  in  the 
downstream  direction.  Consequently,  the  attractor  dimension  grows 
downstream  as  the  instability  amplifies.  Each  downstream  location  in  open 
flow  is  analogous  to  a  different  external  stability  criterion  in  closed  flow. 

Since  each  stability  criterion  in  closed  flow  had  a  characteristic  attractor 
dimension,  so  will  each  downstream  location  in  open  flow  (to  within  the  local 
length  scale)  have  its  characteristic  dimension. 

2.1.3  Romes  to.  Chaos 

There  are  several  well  established  routes  to  chaotic  behavior; 
intermittency,  frequency  locking,  period  doubling  and  the  periodic-quasi- 
periodic-chaotic  sequence.  Each  has  been  observed  in  experiments  conducted 
on  different  types  of  fluid  flows.  Behringer  (1985)  notes  that  the  origins  of 
turbulence  in  convecting  layers  are  usually  due  to  the  nonlinear  interaction 
of  macroscopic  modes  rather  than  microscopic  fluctuations.  He  goes  on  to 
observe  that  "a  strange  attractor  is  a  very  complex  region  of  phase  space,  now 
commonly  associated  with  the  onset  of  turbulence".3 


3Behringer  (1985)  p672. 
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2.1.3. l  Intgrmittgn.cy 


Some  systems  exhibit  a  transition  from  periodic  behavior  (R  <  Rx,  where 
R  is  some  transition  parameter  and  Rx  is  some  critical  value)  to  chaotic 
behavior  (R  >  Rx)  which  is  characterized  by  occasional  bursts  of  noi  e.  For  R 
only  sightly  greater  than  Rx,  there  are  long  intervals  of  periodic  behavior 
interrupted  by  short  bursts  of  noise.  With  increasing  values  of  R,  the  interval 
between  bursts  of  noise  decreases  until  it  eventually  becomes  impossible  to 
distinguish  the  original  underlying  periodic  state.  Behringer  (1985)  notes  that 
intermittency  occurs  when  a  stable  and  unstable  attractor  merge. 

Intcrmittcncy  as  a  route  to  chaotic  behavior  has  been  observed  in 
convection  experiments  (Swinney,  1983)  and  Raylcigh-Benard  convection 
(Behringer.  1985). 

2.1. 3.2  Erequgncy— Locking 

Frequency  locking  is  defined  as  the  transition  from  a  quasi-periodic 
state  to  a  frequency  locked  (periodic)  state  for  some  increasing  value  of  a 
control  parameter.  The  quasi-periodic  state  persists  over  a  wide  range  of  the 
control  parameter,  followed  by  a  well  defined  transition  to  a  chaotic  state. 

Frequency  locking  has  been  observed  as  a  route  to  chaos  in  Rayleigh- 
Bdnard  convection  (Swinney,  1983). 

2. 1.3.3  EgiifldL Doubling 

Period  doubling  occurs  when  a  single  stable  solution  bifurcates  into 
alternating  between  two  stable  solutions  once  a  critical  threshold  is  reached. 
The  solutions  bifurcate  again  as  the  critical  parameter  further  increases.  The 
ratio  between  successive  bifurcations  is  given  by  Feigenbaum's  number. 

Period  doubling  has  been  observed  as  a  route  chaotic  to  behavior  in 
Rayleigh-Benard  convection  (Swinney,  1983;  Behringer,  1985;  Libchaber  et  al., 
1983)  and  shallow  water  waves  (Swinney,  1983). 
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2.1. 3.4  Pcriodic-Quasi-Pgriodic-Chaotic  Sequence 


The  periodic-quasi-periodic-chaotic  sequence  is  defined  as  the 
transition  of  a  system  which  is  periodic  (i.e.,  characterized  by  only  one 
frequency)  to  quasi-periodic  (i.e.,  a  system  characterized  by  two 
incommensurate  frequencies)  to  one  which  is  chaotic  (i.e.,  a  system 
characterized  by  three  or  more  incommensurate  frequencies). 

This  route  to  chaotic  behavior  has  been  observed  in  Couette-Taylor  flow 
(Swinney,  1983). 


2.2  Chaos  in  the  Atmosphere 

There  have  been  a  number  of  studies  of  strange  attractors  in 
meteorological  data.  Most  of  these  have  focused  on  the  longer  time  scales  of  the 
synoptic  and  climactic  range.  Much  of  the  work  on  the  longer  time  scales  has 
been  prompted  by  Lorenz’s  pionecrinj  identification  of  a  strange  attractor  in 
a  model  of  the  general  circulation  (Lorenz,  1963).  Very  little  work  has  been 
done  on  the  time  scales  over  which  waves  decay  to  turbulence. 

2.2.1  Short  Time  Scales 

Tsonis  and  Eisner  (1988)  searched  for  an  attractor  over  very  short  time 
scales  in  vertical  velocity  data  from  the  boundary  layer.  They  estimated  a 
dimension  of  ~  7.3  for  vertical  winds  measured  at  10  m  height  during  the  day 
in  Boulder,  Colorado.  The  data  consisted  of  10  second  averages  of  10  m  vertical 
winds  measured  over  an  11  hour  period  from  1330  -  0030  GMT,  totalling  3960 
points. 

The  estimate  of  the  system  dimension  by  Tsonis  and  Eisner  (1988)  is 
flawed  for  several  different  reasons.  Smith  (1989)  showed  that  over  2.3X101 1 
points  would  be  required  to  accurately  obtain  a  dimension  of  7.3,  while  earlier 
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estimates  of  the  number  of  points  required  to  obtain  the  attractor  dimension 
(i.e.,  10d2-100d2  )  suggest  that  over  107  points  would  be  needed  to  accurately 
obtain  an  attractor  dimension  of  7.3!  It  is  unlikely  that  one  could  specify  the 
dimension  of  a  system  with  eight  degrees  of  freedom  from  an  analysis  of  only 
3960  data  points. 

Tsonis  and  Eisner  (1988)  had  no  initial  hypothesis  as  to  why  there  would 
be  an  attractor  in  the  data  set.  They  speculated  that  the  attractor  was  connected 
with  a  convective  system  described  by  at  least  eight  differential  equations  and 
might  be  related  to  the  Lorenz  system.  In  addition,  this  data  set  was  certainly 
non-stationary,  a  fact  not  considered  by  Tsonis  and  Eisner  (1988).  Any  estimate 
of  an  attractor  dimension  must  consider  stationarity,  else  the  data  suggest  the 
presence  of  a  finite  dimension  strange  attractor  where  there  is  none  (see 
Chapter  III  for  a  more  thorough  discussion). 

Henderson  and  Wells  (1988)  also  used  vertical  velocity  data  from  the 
boundary  layer  to  estimate  the  dimension  of  an  attractor.  Their  data  consisted 
of  the  vertical  velocities  at  10  m  above  the  ground  measured  by  a  sonic 
anemometer  during  the  passage  of  a  thunderstorm  gust  front  over  an  ~  10 
minute  period.  Henderson  and  Wells  (1988)  found  evidence  of  an  attractor  with 
dimensions  between  4.0  and  5.5.  To  accurately  obtain  a  dimension  of  4  or 
greater  would  have  required  more  than  3xl06  points!  While  Henderson  and 
Wells  (1988)  did  not  specify  the  number  of  points  they  used  in  their  analysis, 
the  number  of  points  required  for  an  accurate  dimension  estimate  implies  a 
minimum  sampling  rate  of  5000  Hz.  It  is  unlikely  that  there  were  a  sufficient 
number  of  points  to  determine  the  attractor  dimension. 

2.2.2  Long  Time  Scales 

Most  studies  looking  for  attractors  in  atmospheric  data  consider  longer 
lime  scales.  Fraedrich  (1986)  used  the  Grassberger-Procaccia  algorithm  to 
determine  the  dimension  of  an  attractor  in  surface  pressure  data,  sunshine 
duration  data  and  500  mb  zonal  wave  amplitude  data.  In  a  later  study 
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(Fraedrich,  1987),  he  examined  similar  data  to  determine  the  dimension  of  the 
systems  under  consideration  and  also  investigated  the  predictability. 

The  surface  pressure  data  consisted  of  15  years  (5475  points)  of  surface 
pressure  measurements  made  at  0600  GMT  in  Berlin,  West  Germany.  Fraedrich 

\ 

considered  two  cases:  15  years  of  the  annual  cycle,  and  seasonal  data  sets  made 
over  14  winter  and  15  summer  seasons.  Fraedrich  (1986)  found  no  evidence  of 
an  attractor  with  a  finite  dimension  for  the  continuous  data  set.  In  contrast, 
Fraedrich  (1986)  calculated  dimensions  of  3.2  and  3.9  for  the  winter  and 
summer  seasons,  respectively.  However,  in  another  study  of  the  same  data  set 
published  a  year  later,  Fraedrich  (1987)  found  a  dimension  of  >  6.8-7. 1.  He  gave 
no  reason  for  the  difference  in  the  latter  finding. 

«' 

Fraedrich  (1986)  repeated  these  calculations  for  a  30  year  record  of  the 
number  of  daily  sunshine  hours.  Again,  the  data  set  was  considered  as  two 
separate  cases:  a  30  year  continuous  record  and  separate  winter  (29)  and 
(  summer  seasons  (30).  As  for  the  surface  pressure  data,  Fraedrich  (1986)  found 

no  evidence  for  an  attractor  with  a  finite  dimension  in  the  continuous  record, 
but  estimated  dimensions  of  3.1  and  4.3  for  the  winter  and  summer  seasons, 
respectively. 

Fraedrich  (1986)  repeated  the  analysis  for  10  years  of  500  mb  zonal  wind 
data  at  50°N.  The  continuous  record  did  not  support  evidence  of  a  finite 
dimension  attractor.  Fraedrich  (1986)  calculated  dimensions  of  3  and  3.6  for 

the  winter  and  summer  cases,  respectively. 

i 

On  the  climactic  scale,  Fraedrich  (1987)  used  an  oxygen  isotope  record 
from  deep  sea  core  analysis  to  obtain  a  dimension  of  4.4-4.8  for  an  attractor. 

The  predictability  of  this  attractor  was  between  10000-15000  years. 

(  Grassberger  (1986)  also  searched  for  evidence  of  a  climactic  attractor  in 

oxygen  isotope  ratios  from  deep  sea  cores.  He  found  no  evidence  of  a  finite 
dimensional  attractor  in  the  data  set.  The  small  number  of  data  points  used  in 

\  the  study  prevented  attributing  a  dimension  less  than  10  to  the  system. 
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Xeppenne  and  Nicolis  (1989)  applied  the  Grassberger-Procaccia 
algorithm  to  9000  days  of  500  mb  geopotential  height  records  from  5  stations 
over  western  Europe.  They  calculated  a  dimension  of  7.5  with  a  dispersion  of 
10%  for  the  attractor  in  this  data.  They  discovered  fractal  dimensions  in  the 
data  from  each  station  as  well  as  the  average  of  the  five  stations.  Keppenne 
and  Nicolis  (1989)  used  empirical  orthogonal  functions  to  support  their 
findings. 

These  studies  suffered  from  the  same  shortcomings  that  plagued  the 
shorter  time  scale  investigations.  In  all  cases,  it  is  difficult  to  support  evidence 
of  an  attractor  given  a  limited  number  of  points  in  the  calculation  of  the 
dimension.  In  a  later  section  we  will  describe  the  number  of  points  necessary 
to  accurately  estimate  the  dimension  of  an  attractor. 

Fraedrich's  studies  (Fraedrich,  1986)  of  the  different  synoptic  scale  data 
would  have  required  at  least  74,088  points,  a  number  far  greater  than  the  7300 
actually  used.  The  revised  dimension  for  the  surface  pressure  attractor,  d  >  6.8 
-  7.1  (Fraedrich,  1987),  would  have  required  5,489,031,744  points,  an  even 
larger  figure!  The  dimension  estimate  for  the  climactic  scale  attractor  obtained 
from  the  oxygen  isotope  data  is  similarly  flawed.  To  support  evidence  of  a 
dimension  of  7.5,  Keppenne  and  Nicolis  (1989)  would  have  needed  more  than 
2.3X1011  points  instead  of  the  18,000  they  used.  However,  their  dimension 
estimate  was  supported  by  a  similar  finding  using  a  completely  independent 
technique  and  thus  can  be  given  more  credence. 

Unlike  those  for  short  time  scales,  large  scale  studies  have  a  stronger 
theoretical  basis;  the  work  of  Lorenz  (1963)  shows  the  presence  of  a  strange 
attractor  for  the  synoptic  or  climactic  time  scales.  The  shorter  time  scales  lack 
this  theoretical  underpinning. 
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2.3  Gravity  Waves  and  Gravity  Wave  Breaking 


Internal  atmospheric  gravity  waves  (sometimes  referred  to  as  buoyancy 
waves)  were  first  proposed  by  Hines  (1960)  as  a  mechanism  for  describing 
observations  of  travelling  ionospheric  disturbances  (TID's).  The  gravity  wave 
mechanism  satisfied  some  important  characteristics  of  the  observations: 
upward  propagation  of  the  wave  and  the  increase  of  the  wave  amplitude  with 
height4. 

Gravity  waves  have  frequencies  in  the  range:  f  «  to  «  N,  where  f  is  the 
Coriolis  parameter,  f  =  2Gsin<)>  ~  1.1  x  10-4  s~!  for  the  latitude  of  the  data  (57°  N) 
to  be  used  in  this  study,  and  N  is  the  Brunt-Vaisalia  frequency,  defined  as 


N2  =  |^ 
e3z( 


(2.2) 


where  g  is  the  acceleration  due  to  gravity  and  9  is  the  mean  potential 
temperature.  This  frequency  is  equivalent  to  a  period  of  ~  5  minutes  in  the 
mesosphere  (Andrews  et  al.,  1987). 

There  arc  many  mechanisms  which  lead  to  the  saturation  and 
dissipation  of  atmospheric  gravity  waves.  The  primary  mechanisms  are 
thought  to  be  dynamic  and  convective  instabilities  although  there  are 
competing  mechanisms.  Other  mechanisms  which  limit  wave  growth  are  wave 
dissipation  by  turbulence,  molecular  diffusion,  radiative  damping,  inertial 
instability,  wave  transience  and  the  cascade  of  wave  energy  to  small  scales  via 
nonlinear  wave-wave  interaction  (Fritts  and  Raslogi,  1985).  Fritts  and  Rastogi 


4  Conservation  of  energy  requires  the  amplitude  of  upward  propagating  waves 
to  grow  at  a  rate  proportional  to  [  p(z)]-1/2.  In  the  atmosphere,  density 
decreases  approximately  exponentially  with  height,  i.e.  at  a  rate  proportional 
to  e*z/H,  where  H  is  the  atmospheric  scale  height.  Thus,  upward  propagating 
waves  grow  at  a  rate  of  approximately  ez^2^. 
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(1985)  note  that  wave-wave  interaction  is  the  most  efficient  of  these 
mechanisms. 

In  the  turbulence  theory  of  gravity  wave  saturation,  the  wave 
amplitude  is  limited  by  the  turbulent  eddy  diffusivity  which  originates  from 
the  convective  instability  of  the  wave  itself,  with  convective  overturning 
hypothesized  not  to  occur.  Walterscheid  and  Schubert  (1990)  object  to  this 
theory  of  gravity  wave  saturation  on  two  points.  First,  their  model,  which 
contains  no  assumptions  about  the  eddy  diffusivity,  shows  that  overturning 
does  indeed  occur.  The  upward  propagating  wave  is  not  limited  to  neutral 
stability,  but  instead  develops  highly  unstable  regions  over  certain  phases  of 
the  wave.  The  overturning  causes  localized  convection  which  restores  neutral 
stability.  Second,  the  turbulence  generated  by  the  breakdown  of  the  wave  does 
not  act  to  limit  growth  of  the  wave,  but  is  a  consequence  of  the  nonlinear 
overturning.  As  Walterscheid  and  Schubert  (1990)  point  out,  overturning  and 
wave  saturation  can  occur  even  in  the  absence  of  turbulence.  Upward 
propagating  waves  overturn  when  the  lapse  rate  of  the  mean  potential 
temperature  plus  wave  potential  temperature  becomes  unstable,  i.e., 


a(e  +  e') 

dz 


<0 


(2.3) 


The  unstable  lapse  rate  is  equivalent  to  the  condition  where  the  wave  plus  the 
mean  horizontal  velocity  exceeds  the  phase  speed  of  the  wave,  i.  e.. 


u '  +  u  >  C. 


(2.4) 


Rearranging  this  expression  yields  a  measure  of  the  degree  of  nonlinearity  of 
the  wave,  i.e., 


when  overturning  occurs.  Thus,  gravity  waves  are  highly  nonlinear  when 
saturation  takes  place. 

Dynamical  or  shear  instability  occurs  when  the  wave  amplitude 
becomes  large  enough  that  the  wave  plus  the  mean  velocity  has  a  Richardson5 
(Ri)  number  less  than  1/4.  Fritts  (1982)  showed  that  regions  of  dynamic  and 
convective  instability  are  essentially  the  same,  but  that  convective  instability 
should  occur  first  and  preempt  shear  instability.  Walterscheid  and  Schubert 
(1990)  show  that  convective  instability  does  occur  first  and  preempts  shear 
instability  by  limiting  the  growth  of  the  wave  amplitude  with  height. 

However,  Chimonas  (1986)  showed  that  waves  can  by  dynamically  unstable  for 
any  Richardson  number  if  the  flow  is  tilted.  The  previous  studies  have  all 
assumed  that  the  flow  is  horizontally  stratified. 

Transient  effects  occur  when  wave  breakdown  modifies  the  mean  flow 
and  can  cause  self  acceleration  of  the  wave.  This  can  then  lower  the  height  at 
which  the  wave  breaks.  Walterscheid  and  Schubert  (1990)  believe  wave 
transience  introduces  important  considerations  in  the  wave  breaking  process 
but  that  it  is  not  the  principle  mechanism  behind  gravity  wave  saturation. 


5  The  Richardson  number  is  defined  as: 


Rin 


g  38 
a  Bz 


J $1 


[5%J2  [3%J2 


where  v  is  the  velocity  and  the  other  variables  retain  their  previous 
definitions. 


2.3.1  Conventional  Saturation  Theory 


The  term  gravity  wave  saturation  refers  to  any  process  that  acts  to  limit 
or  maintain  constant  wave  amplitudes  with  altitude.  The  process  occurs  via 
instabilities  or  interactions  arising  from  large  amplitude  wave  motions. 

Gravity  wave  saturation  plays  an  important  role  in  maintaining  the 
mesospheric  circulation  and  temperature  gradient;  gravity  wave  saturation 
provides  the  drag  necessary  to  explain  the  mean  zonal  wind  reversals 
observed  in  the  upper  mesosphere  and  lower  thermosphere  (see  for  example, 
Fritts,  1984;  Holton,  1982,  1983;  Dunl-erton  1982;  Lindzen,  1981). 

Here,  we  provide  a  brief  review  of  linear  saturation  theory  as  proposed 
by  Lindzen  (1981).  Linear  saturation  theory  assumes  that  the  growth  of  the 
amplitude  of  monochromatic  gravity  waves  in  a  horizontally  stratified  flow 
would  be  limited  by  the  appearance  of  convective  instability.  This  would  result 
in  the  production  of  turbulence  and  a  level  of  eddy  diffusion  that  is  just 
sufficient  to  restrain  wave  amplitudes  to  the  unsaturated  limit.  This  theory 

assumes  that  the  gravity  wave  saturation  does  not  affect  wave  propagation  or 

the  wave  characteristics. 

The  basic  equations  in  Cartesian  coordinates  for  an  inviscid  atmosphere 

are: 


pdT  +  Vp'p®  =  0 

dp  p0N2 
-r-'—r—  w  =  0 


dt 


8 


V-v  =  0 


(2.6) 


where  p  is  the  atmospheric  density,  v  is  the  vector  velocity,  w  the  vertical 
component  of  the  velocity,  g  the  acceleration  due  to  gravity,  p  is  the  pressure 
and  N2  is  the  Brunt-VaisSlia  frequency.  The  equations  in  2.6  are  Euler's 
equation,  a  modified  form  of  the  thermodynamic  equation  and  the  continuity 
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equation,  respectively.  The  Boussinesq  approximation  is  implicit  in  this  set  of 
equations,  i.e.,  the  atmosphere  will  be  considered  as  an  incompressible  fluid 
except  in  the  buoyancy  term.  This  approximation  is  not  especially  valid  for  the 
atmosphere  as  a  whole,  but  may  be  valid  for  the  region  near  the  level  of 
gravity  wave  saturation. 

Let  us  apply  the  following  perturbations  to  an  incompressible,  inviscid 
and  adiabatic  atmosphere.  The  background  is  assumed  to  be  hydrostatic  with  a 
zonal  wind  that  varies  with  height,  i.e., 


u  =  uo(z)  +  u'(x,z,t) 
w  =  w'(x,z,t) 
p  =  po(z)  +  p'(x,z,t) 
p  =  po(z)  +  p’(x,z,t)t 

(2.7) 


where  the  primes  indicate  the  perturbation  quantities,  uo  is  the  basic  state 
zonal  wind,  and  po,  po  arc  the  basic  state  density  and  pressure.  The  basic  state 
pressure  and  density  vary  with  height  as 


Po(z)  =  Po(0)e‘z/H 

po(z)  =  p0(0)e'z/H 

(2.8) 


where  po(0)  and  po(0)  are  the  pressure  and  density  at  the  surface  and  H  is  the 
scale  height  of  the  atmosphere. 

Applying  the  perturbations  (2.7)  to  the  set  of  equations  in  (2.6)  and 
neglecting  second  and  higher  order  terms  yields  the  set  of  perturbation 
equations 
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dx 


dw'  ,  dw'  dp'  ■ 
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dp  dp 
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w‘  =  O 


PoN2  . 
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g 
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dx  dz 


We  can  solve  for  w’  first  by  assuming  solutions  of  the  form 


(2.9) 


u’=  u(z)e1(t0l'kx) 
w'  =  w(z)  ei(“t'kx) 
p  =  p(z)  ei((0l_kx) 
p,=  p(z)ei(“‘-kx) 

(2.10) 

and  by  substituting  these  into  the  perturbation  equations  (2.9),  Cancelling  out 
the  exponential  terms  yields  the  set  of  perturbation  equations: 


ipo[co  -  ku0]u  +  po^w  =  ikp 

ipo[o)-ku0]w  =  -^w-pg 

dw  n 
iku  +  -=r-  =  0 
dz 

i[co-ku0]p  -^-w  =  0 


(2.11) 


Solving  this  for  w  yields  the  Taylor-Goldstein  equation  (Fritts,  1984): 


d2w 


N2 


J _ 92U0  ).  2  1  duo  1 


dz2  L(a> - ku0) 2  (®'kuo)  dz2  H(o>-kuo)dz  4H2 


w  =  0 


(2.12) 
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Liiidzen  (1981)  notes  that  for  cases  of  gravity  wave  saturation,  we  can  assume 
that 


1  ^  N2 

4H2  (co-kUo)2 

(2.13) 

which  allows  us  to  drop  the  last  term  on  the  left  hand  side  of  equation  2.12.  The 
d2uo/3z2  and  3uo/9z  terms  can  also  be  dropped  because  the  basic  state  zonal 
wind  is  assumed  to  be  a  slowly  varying  function  of  height.  Scaling  arguments 
also  allow  us  to  neglect  the  k2  term.  This  reduces  equation  2.10  to 


82w 

dz2 


+ 


N2  k2 
(o)-  kuo)2 


w  =  0 


(2.14) 


Near  the  critical  level,  the  denominator,  (to-kuo),  goes  to  zero  because  the 
absolute  value  of  the  horizontal  component  of  the  phase  speed  (to/k) 
approaches  the  background  zonal  wind  speed.  This  zero  creates  a  singularity  at 
the  critical  level,  making  solution  of  equation  2.14  difficult. 

The  WKBJ  approximation  is  used  to  solve  an  equation  of  the  form  of  2.14. 
The  WKBJ  method  is  described  as  follows  (Mathews  and  Walker,  1970).  Given  an 
equation  of  the  form 


82y 

8x2 


+  f(x)y  =  0 


(2.15) 


where  f(x)  is  a  slowly  varying  function  of  x  and  does  not  pass  through  a  zero 
or  other  singularity,  then  solutions  to  equation  2.15  are  of  the  form 


(2.16) 


V7(x)dx 


If  we  define  the  following  quantity, 


?>2,  N2 
(co-kUo)2, 


then  the  solution  to  equation  2.15  becomes 


w  =  A  X 


Xti. 


(2.17) 


(2.18) 


and  is  depicted  in  Figure  2.1. 

Lindzcn  (1981)  uses  this  result  to  show  that  the  condition  for  the 
convective  saturation  of  gravity  waves  is 


fTL 

dz 


(2.19) 


where  T  is  the  perturbation  temperature  and  r  is  the  dry  adiabatic  lapse  rate. 
Other  commonly  used  conditions  for  convective  gravity  wave  saturation  are: 


Po  +  p'l  =  0, 


(2.20) 


Po  +  p'  =0, 


(2.21) 
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(2.22) 


dteo  +  e') 

dz 


=  0 


and  |  u0  +  u  |  =  0. 

(2.23) 

Saturation  conditions  in  2.20  and  2.21  occur  because  you  can  not  have  negative 
pressures  or  densities.  Condition  2.22  is  the  same  as  2.19;  it  is  the  condition  for 
ordinary  convective  instability.  The  last  condition,  2.23,  relates  to  the  phase 
speed  of  the  wave  "catching"  up  to  the  basic  state  wind  speed. 

The  gravity  wave  saturation  model  of  Wtlterschcid  and  Schubert  (1990) 
provides  some  interesting  insights  into  the  mechanism  behind  gravity  wave 
saturation.  Two  dimensional  and  fully  nonlinear,  the  model  makes  no 
significant  assumptions  about  the  gravity  wave  saturation  mechanism;  the 
atmosphere  in  the  model  is  compressible  and  non-hydrostatic.  Most 
importantly,  the  model  indicates  that  nonlinear  growth  of  the  gravity  wave 
creates  regions  of  overturning  prior  to  saturation.  The  large  unstable 
potential  temperature  gradients  allow  the  development  of  small  scale  cellular 
convection  which  causes  the  breakdown  of  the  gravity  wave.  Turbulence  is  an 
end  product  of  the  wave  breakdown  via  the  decay  of  the  cellular  convection 
rather  than  the  cause  of  the  wave  breakdown.  Walterscheid  and  Schubert 
(1990)  cite  the  laboratory  work  of  Delisi  and  Corcos  (1973)  as  support  for  this 
conclusion  from  their  model. 

2.3.2  Slantwise  Static  Instability  Theory 

Hines  (1988)  objects  to  the  linear  saturation  mechanism  proposed  by 
Lindzen  (1981)  on  two  points:  the  spectrum  of  waves  and  the  vertical 
gradients.  In  raising  these  objections,  Hines  suggests  that  slantwise  static 
instability  is  a  less  demanding  mechanism  for  gravity  wave  breaking. 

Waves  in  the  middle  atmosphere  are  not  represented  by  a  single 
dominant  wave  number  but  instead  by  a  spectrum  of  wave  numbers.  Hines 
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(1988)  notes  that  this  does  not  affect  the  linear  saturation  theory  of  wave 
breaking  per  se,  but  it  "does  open  the  way  to  patchiness  rather  than  a  laminar 
deposition  of  turbulence"6.  The  instability  criterion  for  a  single  wave  mode, 
i.e.,  where  the  phase  speed  of  the  perturbation  equals  that  of  the  background 
flow,  does  not  impose  nonlinearity  on  the  system.  It  would  impose  large 
nonlinearities  when  two  waves  have  similar  amplitudes  but  different  wave 
vectors.  Hines  (1988)  states  that  the  nonlinear  interaction  may  leach  away 
wave  energy  before  the  wave  becomes  unstable. 

Hines  (1988)  also  objects  to  the  conventional  saturation  theory  because 
it  assumes  the  background  flow  is  strictly  horizontally  stratified  and  only 
considers  vertical  gradients  in  the  evaluation  of  the  stability.  Hines  notes  that 
this  restriction  is  more  mathematical  in  nature  but  is  not  justified  physically. 
Hines  (1988)  examines  the  stability  criteria  for  gradients  which  are  no  longer 
strictly  vertical. 

Hines  (1988)  notes  that  the  criterion  for  the  onset  of  instability  in  a 
horizontally  stratified  atmosphere  is  given  by 


>o 

0  9z 


(2.24) 


where  g  is  the  acceleration  due  to  gravity,  8  is  the  mean  potential  temperature. 
This  is  the  negative  of  the  Brunt-VaisSlla  frequency  wb2  (defined  as  N2  in  the 
previous  section)  and  as  such  would  appear  in  the  form 


gicobt 


(2.25) 


6  Hines  (1988),  p!269. 


in  the  solution  of  the  relevant  linearized  equations.  If  equation  2.25  were 
negative,  i.e.,  unstable  conditions,  then  the  solutions  would  take  the  form 


e 


(2.26) 


where  tz  is  the  e-folding  time  for  the  growth  of  any  instabilities.  This 
mechanism  will  be  referred  to  hereafter  as  vertical  static  instability.  If  the 
atmosphere  is  not  horizontally  stratified,  i.c.,  the  gradient  of  the  potential 
temperature  is  no  longer  constrained  to  the  vertical,  the  new  criterion  for 
instability  is  given  by 
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cos  £  >  0 


(2.27) 


where  £  is  the  angle  off  the  vertical  for  any  parcel  motion  oriented  along  the  s 
axis  as  shown  in  Figure  2.2.  Any  parcel  motion  confined  to  the  shaded  region 
in  the  figure  would  be  unstabie.  By  analogy,  the  c-folding  time  for  growth  of 
an  instability  for  this  type  of  stratified  atmosphere  would  be  given  by 


e 


(2.28) 


The  e-folding  time  for  growth  of  the  instability  may  be  long  for  certain  values 
of  C- 

Hines  (1988)  concludes  that  turbulence  is  far  more  likely  to  develop 
from  slantwise  static  instability  than  for  vertical  static  instability,  even 
though  it  may  be  a  slower  mechanism.  Since  gravity  waves  produce  potential 
temperature  gradients  that  are  not  necessarily  constrained  to  the  vertical 
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plane,  this  mode  of  instability  is  the  more  likely  to  generate  turbulence.  Under 
most  conditions,  the  normalized  wave  amplitude  necessary  for  the 
development  of  turbulence  from  slantwise  static  instability  is  «1;  the 
normalized  wave  amplitude  necessary  for  the  development  of  turbulence  for 
vertical  static  instability  is  defined  as  1. 

Hines  also  concludes  that  the  turbulence  spectrum  from  slantwise  static 
instability  would  be  highly  anisotropic,  with  much  stronger  horizontal 
motions  than  vertical  motions.  This  might  approach  two  dimensional 
turbulence  in  the  limit.  The  production  of  turbulence  from  slantwise  static 
instability  is  more  likely  for  "short"  vertical  wavelengths  than  from  vertical 
static  instability.  Hines  (1988)  notes  that  for  vertical  wavelengths  of  6  km  and 

a  buoyancy  period  of  5  minutes,  vertical  static  instability  requires  vertical 
velocities  of  20  m  s*1  for  saturation  to  occur.  This  is  not  consistent  with 
observations.  However,  because  slantwise  static  instability  requires  much 
smaller  wave  amplitudes  for  saturation,  this  mechanism  can  produce 

saturation  consistent  with  observations  of  vertical  wind  velocities  in  the 
mesosphere. 

The  one  difficulty  with  finding  a  strange  attractor  associated  with  the 
slantwise  static  instability  mechanism  is  that  the  smaller  wave  amplitudes 
which  produce  saturation  and  turbulence  may  allow  the  gravity  wave  to 

remain  fairly  linear.  The  analysis  by  Hines  (1988)  was  done  for  a 
monochromatic  wave  to  simplify  the  mathematics  and  is  limited  by  this 
assumption.  Nonlinear  wave-wave  interaction  will  affect  both  saturation 
mechanisms  equally.  Without  nonlinear  wave-wave  interaction,  the  slantwise 

static  instability  mechanism  reduces  the  probability  of  finding  an  attractor  in 
the  saturation  of  gravity  waves.  However,  the  "patchiness"  of  the  turbulence 
created  by  saturation  of  gravity  waves  with  a  spectrum  of  wave  numbers 
opens  the  way  for  the  investigation  of  the  fractal  structure  of  turbulence. 
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2.3.3  Nonlinear  Mechanisms 

There  are  three  "resonant  triad"  interactions  which  act  on  gravity 
waves:  elastic  scattering,  induced  diffusion  and  parametric  subharmonic 
instability  (Fritts  and  Rastogi,  1985).  Of  these  mechanisms,  Fritts  and  Rastogi 
feel  parametric  subharmonic  instability  acts  most  efficiently  in  transferring 
wave  energy  between  waves  of  very  different  scales.  Parametric  subharmonic 
instability  transfers  energy  from  large  scale  waves  to  two  small  scale  waves  at 
half  the  frequency  of  the  larger  wave.  Elastic  scattering  converts  the  incident 
wave  into  a  reflected  wave  by  scattering  off  a  vertical  shear  in  the  wind  and 
tends  to  make  the  vertical  wave  spectrum  symmetric.  Elastic  scattering  can  be 
rapid,  depending  on  the  vertical  wave  number  (Ych  and  Liu,  1985).  In  induced 
diffusion,  the  wave  action  density  diffuses  in  wave  space  if  two  waves  with 
nearly  identical  wave  vectors  interact  with  the  vertical  shear  of  a  smaller 
wave.  Ibrahim  (1987)  suggests  that  induced  diffusion  can  be  sufficient  for 
gravity  wave  saturation. 

Fritts  and  Rastogi  (1985)  note  that  parametric  subharmonic  instability 
acts  on  waves  of  all  amplitudes  and  may  exchange  energy  among  waves  at 
amplitudes  less  than  that  required  for  convective  or  dynamic  instability. 
However,  it  works  best  at  high  vertical  wave  numbers  and  small  intrinsic 
frequencies.  Fritts  and  Rastogi  do  not  believe  that  parametric  subharmonic 
instability  competes  effectively  with  convective  or  dynamic  instability  among 
higher  frequency,  larger  scale  gravity  wave  motions  as  a  mechanism 
explaining  gravity  wave  saturation. 

Parametric  subharmonic  instability  transfers  the  energy  in  moderate- 
to-large  scale  waves  to  small  scale  waves  at  one  half  the  frequency.  Yeh  and 
Liu  (1985)  note  that  conservation  of  wave  number  and  conservation  of 
frequency  restrict  this  mechanism  to  those  waves  with  an  elevation  angle  of 
60°  or  greater.  The  time  scale  for  parametric  subharmonic  instability  varies 
with  the  inverse  square  of  the  vertical  wavenumber. 
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Yeh  arid  Liu  (1985)  state  that  if  wave  growth  occurs  over  an  amount  of 
time  equal  to  a  few  periods  or  less,  the  rionlinear  interactions  are  no  longer 
considered  "weak"  and  other  processes,  such  as  convective  or  dynamic 
instability,  may  become  important.  Furthermore,  the  interactions  in  their 
study  only  considered  the  energy  transferred  in  the  vertical  plane  and 
neglected  energy  transferred  in  the  horizontal  plane. 

Nonlinear  wave-wave  interaction  can  act  to  limit  the  growth  of  gravity 
waves  by  transfering  energy  from  larger  amplitude  waves  to  smaller  waves. 
This  mechanism  will  not  necessarily  lead  to  gravity  wave  breaking.  There  may 
be  no  transition  from  a  laminar  state  to  a  turbulent  one  and  therefore  no 
potential  strange  attractor  in  the  flow.  There  may  be  a  strange  attractor 
associated  with  nonlinear  wave-wave  interaction,  but  it  may  be  a  function  of 
the  underlying  spectrum  of  waves  instead  of  a  function  of  the  transition  from 
laminar  to  turbulent  flow.  The  latter  is  more  likely  to  be  a  universal  and 
repeatable  behavior,  whereas  the  former  will  only  describe  a  particular 
packet  of  gravity  waves  and  may  not  be  found  again. 
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Figure  2.1  Nonvertical  orientation  for  the  potential  tempctature  gradient. 

Interchanges  of  air  parcels  along  axes  within  the  shaded  regions 
such  as  the  s  axis  are  unstable  (adapted  from  Hines,  1988). 
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CHAPTER  III 


THE  METHODS  OF  ANALYSIS 


This  chapter  describes  the  methods  that  will  be  used  to  analyze  the  data. 
The  bulk  of  the  chapter  is  devoted  to  a  description  and  derivation  of  the 
Grassberger-Procaccia  correlation  integral  algorithm,  since  this  technique  is 
not  commonly  used  on  meteorological  data.  A  discussion  of  the  strengths, 
weaknesses  and  difficulties  that  arise  in  employing  this  algorithm  follows  the 
derivation.  The  remainder  of  the  chapter  provides  a  brief  overview  of  the 
more  conventional  analysis  tools  (autocorrelation  and  power  spectrum)  that 
will  be  used  to  support  the  analysis  of  the  system  dimension. 


3-1  The  Grassberger-Procaccia  Correlation  Integral  Algorithm 

There  arc  a  number  of  techniques  described  in  the  literature  for 
estimating  the  dimension  of  a  strange  attractor.  Box  counting  algorithms  yield 
an  estimate  of  the  capacity  dimension,  commonly  referred  to  as  the  fractal 
dimension  (see,  for  example,  Liebovitch  and  Tibor,  1989).  The  nearest 
neighbor  method  developed  by  Badii  and  Politi  (1987)  is  another  approach  to 
estimating  the  dimension  of  a  system.  Yet  another  dimension  estimate  can  be 
obtained  from  singular  systems  analysis  (see  Broomhead  and  King,  1986; 

Albano  et  al.,  1988;  Vautard  and  Ghil,  1989).  The  spectrum  of  Lyapunov  (or 
characteristic)  exponents  for  an  attractor  can  be  calculated  and  related  to  the 
dimension  of  the  system  (see  Packard  et  al.,  1980;  Froehling  et  al.,  1981;  Roux  et 
al.,  1983;  Wolf  et  al.,  1985).  However,  the  most  frequently  used  method  to 
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calculate  the  dimension  of  a  system  is  the  correlation  integral  algorithm 
(Grassbergcr  and  Procaccia,  1983;  Grassberger  and  Procaccia,  1984). 

Each  technique  has  its  strengths  and  weaknesses.  The  box  counting 
technique  is  computationally  intensive  although  Liebovitch  and  Tibor  (1989) 

developed  a  faster  and  more  efficient  algorithm.  The  singular  systems 
approach  is  difficult  to  implement  although  it  is  often  used  as  a  check  on  the 
other  methods  (Broomhead  and  King,  1986;  Albano  et  al.,  1988;  Vautard  and 
Ghil,  1989).  The  nearest  neighbor  approach  is  relatively  new  and  has  not  been 
widely  used;  its  strengths  and  weaknesses  have  yet  to  be  thoroughly  examined 
in  the  literature.  Calculating  the  spectrum  of  Lyapunov  exponents  often 
presupposes  some  knowledge  of  the  attractor.  The  Grassberger-Procaccia 
algorithm  has  the  advantage  of  being  easily  implemented  and  calculated. 

Because  of  its  frequent  use,  the  limitations  of  the  Grassberger-Procaccia 
algorithm  have  been  widely  investigated  and  described  in  the  literature. 

The  Grassberger-Procaccia  algorithm  (along  with  subsequently 
developed  variations)  will  be  used  in  this  study.  It  is  easy  to  understand  and 
program.  It  makes  no  assumptions  about  the  presence  of  a  strange  attractor  in 
the  signal  and  requires  no  a  priori  knowledge  of  its  nature  or  structure. 
Somewhat  computationally  intensive  (the  computer  time  increases  as  the 
square  of  the  number  of  points  in  the  data  set),  the  Grassberger-Procaccia 
algorithm  is  less  demanding  than  the  box  counting  algorithm.  Despite  these 

advantages,  this  algorithm  has  several  drawbacks  which  will  be  discussed  in 
detail  later  in  the  chapter. 

Before  deriving  the  Grassberger-Procaccia  algorithm  and  discussing 
some  of  its  advantages  and  disadvantages,  let  us  define  the  different  concepts 
of  dimension  and  the  method  for  building  a  phase  space  portrait  from  a  single 
data  set. 
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3.1.1  Concepts  of  Dimension 


Three  different  dimensions  are  often  used  to  describe  a  set1:  the 
capacity  dimension,  the  Hausdorff  dimension  and  the  information  dimension. 
These  are  summarized  in  Table  3.1  (adapted  from  Fanner  et  al.,  1983).  The 
Hausdorff  and  capacity  dimension  are  metric  dimensions,  a  concept  of 
dimension  on  which  a  sense  of  distance  is  defined.  The  information  dimension 
is  a  probabilistic  dimension  based  on  the  natural  measure  -  the  relative 
probability  of  different  regions  of  the  attractor  as  obtained  from  time 
averages  (Farmer  et  al.,  1983). 


Table  3.1  Dimension  definitions. 


Name 

Symbol 

Generic  name 

Capacity  dimension 

dc. 

fractal 

Hausdorff  dimension 

dh 

Information  dimension 

di 

The  capacity  dimension  is  defined  as 


dc  = 


Hm  lQgN(£) 
s-0  log(  l/e) 


(3.1) 


where  N(e)  is  the  number  of  cubes  with  sides  of  length  e  needed  to  cover  the 
set  of  points  (Barnsley,  1988).  This  definition  of  dimension  is  the  basis  of  the 
box  counting  algorithm. 

The  Hausdorff  dimension  is  more  complicated  than  the  capacity 
dimension..  Its  definition  is  similar  to  that  of  the  capacity  dimension  but  the 
cubes  used  to  cover  the  set  can  be  of  variable  length.  We  will  not  give  a  more 

1  A  set  can  be  a  collection  of  points,  a  geometric  object  or  a  time  series  of  data. 
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precise  definition  but  note  instead  that  for  most  attractors  the  Hausdorff 
dimension  is  equal  to  the  capacity  dimension  (Farmer  et  al.,  1983). 

The  information  dimension  is  a  generalization  of  the  capacity 
dimension;  it  takes  into  account  the  relative  probability  of  the  cubes  used  to 
cover  the  set.  It  is  defined  as 


,  ,•  1(e) 

di=  hm  — rrn 
e-*o  iog(  V e) 


where  1(e)  is  the  information  for  length  scale  e  defined  as 


(3.2) 


N(e) 

I(e)  =  -X  Pi  log  Pi 

i=l 

(3.3) 

and  pi  is  the  probability  that  the  attractor  trajectory  visits  the  ith  cube.  If  all 
the  cubes  are  visited  with  equal  frequency  then  the  probability  that  the  ith 
cube  is  visited  is  given  by 


Pi  = 


_1_ 

N(e) 


and  therefore  the  information  can  be  written  as 


(3.4) 


1(e)  =  log  N(e), 

(3.5) 

Thus,  the  capacity  dimension  is  equal  to  the  information  dimension  for  a 
completely  homogeneous  attractor.  In  general,  attractors  are  not 
homogeneous  and  the  cubes  arc  not  visited  with  equal  frequency.  Thus,  for  an 
inhomogeneous  attractor, 
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1(e)  <  log  N(e) 


(3.6) 


and  consequently  the  capacity  dimension  is  always  greater  or  equal  to  the 
information  dimension,  i.e., 


dc>  di 


(3.7) 


Information  theory  gives  a  very  specific  meaning  to  1(e).  It  is  the 
amount  of  information  necessary  to  specify  a  system  to  within  accuracy  e  (e  > 

0).  Alternately,  it  can  be  thought  of  as  the  amount  of  information  obtained  by 
making  a  new  measurement  with  an  uncertainty  e. 

1.1.2  Building  Phase  Space  Vectors 

Originally,  there  was  no  universally  accepted  method  for  constructing 
phase  space  vectors  from  a  time  series.  Packard  et  al„  (1980)  pointed  out, 

"...there  is  no  universally  applicable  method  of  phase  space  construction, 

though  the  natuie  of  the  phenomenon  might  suggest  possible  alternatives."2 

However,  Brandstater  et  a!.,  (1983)  suggested  that  phase  space  portraits  can  be 
constructed  by  lagging  the  original  time  series  by  an  arbitrary  amount  to 

obtain  the  second  phase  space  dimension.  This  technique  of  constructing  a 
phase  space  representation  of  the  data  is  often  referred  to  as  Taken's  method  of 
delays  and  is  now  almost  universally  used. 

Fraser  and  Swinney  (1986)  provide  a  clear  explanation  of  building  a 
multi-dimensional  phase  portrait  from  a  single  time  series.  A  scalar  time  series 
s(t)  can  be  expanded  into  a  m  dimensional  phase  space  vector  x(t)  by  using 
time  delays  t,  as  follows 


2  Packard  et  al.  (1980),  p713. 
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X  ( t)  =  {  X0(t),  X  l(t),  .  Xm.i(t)}  > 


(3.8) 


where 


xn(t)  =  s(t  + n  x);  n=0,l ,2,  ....  m-1 

Here  s  is  the  original  time  series,  x  the  specified  time  lag  and  xo,  xj,  ....  xm-i  are 
the  individual  components  of  an  m  dimensional  phase  space  vector  x.  Fraser 
and  Swinncy  (1986)  note  that  for  an  infinite  amount  of  noise  free  data,  the 
time  delay  x  can  be  arbitrary.  However,  for  noisy  or  limited  data,  a  small  time 
delay  x  may  make  the  components  Xo(t)  and  X](t)  indistinguishable  and  all 
trajectories  appear  to  be  on  a  line  xo=xi.  To  avoid  this  problem,  the  time  delay  x 
must  be  chosen  to  make  the  vectors  xo(t)and  Xj(t)  as  independent  as  possible. 
Techniques  in  determining  the  proper  choice  for  the  time  delay  will  be 
discussed  in  later  sections. 

3.1.3  Trajectories. -in  Phase  Space 

The  trajectory  in  phase  space  is  said  to  follow  an  attractor  if  its  orbits 
rapidly  return  to  this  subset  (i.e.  the  attractor)  after  finite  perturbations 
(Swinncy,  1983).  Large  perturbations  could  send  the  orbit  out  of  the  basin  of 
attraction.  An  attractor  is  labelled  strange  if  nearby  orbits  diverge 
exponentially  on  average.  This  condition  is  sometimes  referred  to  as  "sensitive 
dependence  on  initial  conditions"  (Swinney,  1983). 

3.1.4  Derivation  of  the  Grassbereer-Procaccia  Algorithm 

Grassberger  and  Procaccia  (1984)  note  that  "two  of  the  most  basic 
properties  of  dissipative  chaotic  systems  are  related  to  information:  the 
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Kolmogorov  (or  'metric')  entropy  K  and  the  Renyi-Balatoni  information 
dimension  a."3  Both  K  and  a  relate  to  the  information  I(e,  T)  gained  by 
observing  the  trajectory  of  a  system  with  precision  e  during  a  finite  amount  of 
time  T.  Grassberger  and  Procaccia  (1984)  define  the  precision  e  as  the 
uncertainty  in  the  measurements  of  any  of  the  coordinates  of  the  vector  x. 

The  Kolmogorov  entropy  is  defined  as 


K  =  lim  lim 

(3.9) 

The  definition  of  the  Kolmogorov  entropy  requires  making  a  very  long 
series  of  observations  as  can  be  s,en  from  the  limit  placed  on  time. 

Grassberger  and  Procaccia  (1984)  show  that  since  the  time  limit  is  taken  first, 
equation  3.9  implies  that  the  information  for  a  given  precision  e  increases 
linearly  with  time  and  that  the  rate  of  increase  tends  towards  a  finite  constant 
for  infinite  precision  (i.c.,  infinitely  small  error). 

Furthermore,  Grassberger  and  Procaccia  (1984)  compare  this  to  an 
ordered  system  where 


lim  =  0 

T->~  T 


and  systems  with  random  noise  where 


(3.10) 


lim  x  In  ( 1  /  e) 

T  — >  oo  1 


— > 
E— >0 


oo 


(3.11) 


3  Grassberger  and  Procaccia  (1984),  p35.  Note  that  o  is  the  same  as  the 
previously  defined  information  dimension  di. 
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This  leaves  us  with  the  following  definitions  for  the  behavior  of  the 
Kolmogorov  entropy: 

1.  K  =  0  for  ordered  systems, 

2.  K  =  °°  for  random  systems,  and 

3.  K  =  a  finite  constant  for  systems  characterized  by  a  strange  attractor. 

To  extend  this  concept  to  a  time  series  of  data,  Grassberger  and  Procaccia 
(1983)  introduce  a  new  quantity  K2  which  has  the  following  properties: 

1.  K2>0, 

2.  K2<K, 

3.  K2  =  00  for  random  systems,  and 

4.  K2  *  0  for  chaotic  systems. 

Grassberger  and  Procaccia  (1983)  note  that  K2  >  0  is  a  sufficient  condition  for 
chaos.  The  quantity  K2  can  be  calculated  in  the  following  manner. 

Grassberger  and  Procaccia  (1983)  define  a  new  quantity,  C(e),  which  is  total 
probability  that  a  random  pair  of  points  on  the  attractor  will  fall  into  the  same 
cube  of  size  e  in  phase  space.  Grassberger  and  Procaccia  (1984)  note  that  this 
probability  scales  as 


C(e)  =  ev 

e— >0 


(3.12) 


where  v  is  called  the  correlation  exponent.  Furthermore,  v  also  approximates 
the  fractal  dimension,  dc,  of  the  attractor  (Grassberger  and  Procaccia,  1983). 

For  a  time  series  of  data  \ Xj /j=i  .where  Xj  =  X(t  =  ix),  the  correlation 
integral,  Cm(e),  can  be  calculated  from  the  following  (Grassberger  and 
Procaccia,  1983): 
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Cm(e)=  lim  ■ 

N — >oo 


N2 


x  {number  of  pairs  of  points  (n,k)  with  Xn-Xk  <£ 


(3.13) 


Here  N  is  the  total  number  of  points  in  the  time  series  and  the  time  delay  t 
equals  XAt  for  some  lag  X.  The  subscript  m  is  the  embedding  dimension  and 
must  be  greater  than  or  equal  to  F,  the  number  of  degrees  of  freedom  of  the 
attractor.  The  Whitney  embedding  theorem  states  that  it  is  possible  to  embed  an 
m  dimensional  geometric  object  arbitrarily  in  a  2m+l  dimensional  space. 

Equation  3.13  may  be  rewritten  more  precisely  by  replacing  the  norm 
with  an  explicit  expression  for  the  Euclidean  norm  in  a  m  dimensional  phase 
space: 


:m(e)=  lim 

N — >oo 


N2  \ 


x  !  number  of  pairs  (n,k)  with 


m- 1 


I  |xn+j  -  X|(+j 


L  i=0 


(3.14) 


This  should  be  the  same  as  (Grassberger  and  Procaccia,  1983) 


Cm(e)  *  £vexp(-mxK2) 

m  — 
e— *0 

(3.15) 

This  now  gives  us  a  simple  way  to  calculate  the  fractal  dimension  as  well 
as  estimating  the  lower  bound  on  the  Kolmogorov  entropy..  If  we  plot  the 
natural  logarithm  of  Cm(e)  as  a  function  of  the  natural  logarithm  of  e  for 
increasing  values  of  embedding  dimension  m.  we  should  get  a  series  of 
straight  lines  whose  slope  is  v,  the  fractal  dimension  of  the  attractor  and 
which  are  displaced  from  one  another  by  the  factor  -nvrK2-  An  example  of  the 
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correlation  integral  is  shown  in  Figure  3.1.  Calculating  a  value  for  K2  will 
quickly  tell  us  whether  the  data  is  ordered,  chaotic  or  random. 

In  practice,  the  region  of  constant  slope  will  only  be  valid  over  a  limited 
range  of  e  which  is  often  called  the  "scaling"  region.  This  gives  equation  3.15 
the  form 

C(e)  =  x(e)ev 

(3.16) 

where  x(e)  is  a  possibly  oscillatory  function  of  0(1)  (Smith,  1988).  The  structure 
of  x(£)  is  generated  by  the  sparse  or  empty  regions  (lacunae)  of  the  set.  Some 
of  the  oscillations  in  x(e)  are  also  generated  from  noise  and  fluctuations  due  to 
a  finite  number  of  points. 

The  Kolmogorov  entropy  of  the  attractor  can  be  approximated  by 
examining  the  following 


K2'™(E)  =  >lC„.,(e) 


Cm  (£) 


(3.17) 


and  then 


lim  K2,m(e)  =  K2 

m— >00 
E— >0 


(3.18) 


The  Kolmogorov  entropy,  along  with  the  correlation  dimension,  can  be  used  to 
test  for  strange  attractors  in  a  set  of  data.  An  example  of  using  the 
Grassberger-Procaccia  approximation  of  the  Kolmogorov  entropy  to  determine 
the  type  of  behavior  of  a  system  is  shown  in  Figure  3.2. 

Grassbcrger  and  Procaccia  (1983)  demonstrated  this  technique  on  the 
Mackey-Glass  delay  differential  equation  and  th<  Henon  attractor  and  found 


that  this  technique  works  for  fractal  dimensions  as  large  as  7.5  with  ~  30000 
data  points.  However,  later  research  with  the  correlation  integral  suggests 
that  its  only  practical  for  determining  the  dimension  of  systems  with 
dimensions  on  the  order  of  4  or  less  (Smith,  1988). 

3.1.5  Generalized  Correlation  Integral 

The  technique  developed  by  Grassberger  and  Procaccia  (1983;  1984)  has 
been  extended  to  determine  the  generalized  entropy  and  generalized 
dimension  (Hcntschcl  ct  al.,  1983;  Pawelzik  and  Schuster,  1987;  Grassberger, 
1985).  Such  generalized  quantities  reveal  important  information  on  the 
structure  of  the  attractor. 

Atmanspacher  et  al.,  (1988)  notes  that  the  quantity  for  characterizing 
an  attractor  as  a  metric  structure  is  its  dimension.  Traditionally,  the  concept  of 
dimension  has  been  limited  to  purely  integer  values.  However,  the  dimension 
can  take  on  non-integer  values  for  chaotic  (or  strange)  attractors.  Attractors 
with  purely  integer  dimensions  correspond  to  regular  (i.e.  stationary,  and/or 
periodic)  processes.  The  concept  of  a  fractal  (i.e.  non-integer)  dimension, 
d  <  m,  of  an  attractor  in  a  m  dimension  phase  space  can  be  derived  from 
information  theory.  The  information  dimension,  dj,  describes  how  information 
1(e)  scales  with  varying  spatial  resolution  e  as  previously  defined  in  equation 
3.2. 

One  way  to  obtain  the  information  I  is  to  break  the  attractor  up  into  m 
boxes  of  size  e.  The  probability  that  a  point  on  the  attractor  falls  into  the  ith 
box  is  given  by 


Pi  = 


Ni 

N 


(3.19) 
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where  Ni  is  the  number  of  points  in  the  ith  box  and  N  is  the  total  number  of 
points  on  the  attractor. 

Pawclzik  and  Schuster  (1987)  define  the  generalized  information  of 
order  q  as 


I<q)  = 


(3.20) 


A  continuous  spectrum  of  dimensions  of  order  q  can  be  defined  by  substituting 
into  the  original  definition  of  the  information  dimension  (equation  3.2), 


d(q> 


log  £  p? 

=  1  i  =  1 
q-1  log(e) 


(3.21) 


Some  of  the  most  frequently  encountered  dimensions  are:  dW,  the  Hausdorff 
dimension  (previously  referred  to  as  dh  );  d^1),  the  information  dimension 
(previously  referred  to  as  di);  and  d^2),  the  correlation  dimension 
(Atmanspacher  et  al.,  1988).  Furthermore,  Atmanspacher  et  al..  (1988)  observe 
that 


d(q)<  d(ci  1  if  q'  <  q 

(3.22) 

The  equality  holds  only  for  completely  homogeneous  probability  distributions, 
i.e.,  Pi  =  1  /  N.  The  more  an  attractor  “bunches”  up  (i.e.  spends  more  time 
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visiting  a  particular  region  of  phase  space),  the  less  homogeneous  the 
probability  distribution.  Differences  arise  between  dimensions  of  different 
order  q  because  of  the  degree  of  inhomogeneity  of  the  attractor,  i.e.,  the 
degree  to  which  the  boxes  are  visited  with  unequal  frequency. 

The  correlation  integral  method  proposed  by  Grassberger  and  Procaccia 
(1983)  is  based  on  the  correlation  dimension  v.  This  has  been  extended  to  a 
dimension  of  arbitrary  order  q  as  follows 


where  H  is  the  Heaviside  step  function  (H(x)=0  if  x<0,  H(x)  =  1  if  x>0),  and  e  is 
the  size  of  the  box  (Hentschcl  et  al.,  1983;  Pawelzik  and  Schuster,  1987; 
Grassberger.  1985).  Notice  that  this  reduces  to  the  original  expression  of 
Grassberger  and  Procaccia  (1983)  for  order  q=2.  By  analogy,  the  generalized 
entropy  of  order  q  may  now  be  written  as 


K(  h )  =  lim  lim  (=J-lnCt9)(£,N)) 

E— >0  n— »°°  ^ 

(3.24) 

These  results  can  be  used  for  a  series  of  single  observations  evenly  spaced  in 
time  by  use  of  Taken's  method  of  delays.  This  lets  us  write  the  generalized 
correlation  integral  as 


C(mq)(e,  N)=  lim 

‘ 5  E— >o 


m  - 1 

^  (  xi  +  k 
k  =  0 


Xj  +  k) 


-  i 

“1 

q-i" 

2 

J 

_ 

_ 

1 

q-1 


(3.25) 
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Pawelzik  and  Schuster  (1987)  point  out  that  this  technique  is  only 
slightly  more  computationally  intensive  than  the  original  correlation 
algorithm  proposed  by  Grassberger  and  Procaccia  (1983).  Most  of  the 
computational  effort  is  spent  in  counting  the  number  of  the  pairs  of  points; 
raising  the  interior  sums  to  the  various  powers  represents  only  a  small 

additional  burden. 

Once  a  spectrum  of  generalized  entropies  have  been  generated  for  a 
time  series,  we  can  determine  the  spectrum  of  dynamical  fluctuations  around 
the  Kolmogorov  entropy.  This  spectrum  can  be  used  to  deduce  properties  and 

structure  of  the  attractor  (Atmanspacher  et  al.,  1988)  or  applied  to  a 

description  of  turbulence  (Chhabra  et  al.,  1989;  Meneveau  and  Nelkin,  1989). 


3.2  Limitations  of  the  Grasshergcr-Proeaccia  Algorithm 

The  Grassberger-Procaccia  algorithm  docs  have  a  number  of 
weaknesses  which  must  be  addressed.  First,  noise  in  the  signal  can  yield 
misleading  estimates  of  the  dimension.  Second,  an  inappropriate  choice  of  the 
time  delay  in  constructing  the  phase  space  vectors  can  also  yield  incorrect 
estimates  of  the  attractor  dimension  and  may  suggest  the  presence  of  a  chaotic 
attractor  where  there  is  not  one.  Furthermore,  a  limited  or  non-stationary  data 
set  can  introduce  errors  and  require  the  use  of  different  norms  in  calculating 
the  distances  between  pairs  of  points. 

3.2.1  Noise 

Noise  affects  length  scales  over  a  range  on  the  order  of  magnitude  of 
the  standard  deviation  of  the  noise  (Ben-Mizrachi  et  al.,  1984).  This  leads  to  a 
noise  length  scale  region  where  noise  scales  as  the  embedding  dimension 
(Ben-Mizrachi  et  al.,  1984;  Theiler,  1987;  Franaszek,  1987).  Thus,  noise  is 
proportional  to  £m  and  has  a  slope  on  the  In  C(e,N)  vs  In  e  plots  that  equals  m, 
the  embedding  dimension.  The  presence  of  noise  reduces  the  scaling  region  of 
other  signals;  in  systems  with  a  low  signal  to  noise  ratio,  the  scaling  region 
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may  even  disappear.  The  magnitude  or  amount  of  noise  in  the  signal  can  be 
determined  from  the  break  in  the  slopes  in  the  ln-ln  plots  of  the  correlation 
integral;  this  break  is  often  referred  to  as  a  "knee".  The  presence  of  noise  in 
the  signal  violates  the  limit  in  equation  3.12;  you  can  not  take  the  limit  as  e 
goes  to  zero  because  the  attractor  is  not  clearly  defined  for  length  scales  on 
the  order  of  the  magnitude  of  the  noise. 

Furthermore,  a  limited  number  of  data  points  has  much  the  same  effect 
as  noise  on  the  slopes  from  the  correlation  integral.  This  occurs  because 
random  noise  never  completely  saturates  an  infinite  dimensional  phase  space; 
a  limited  data  set  may  not  cover  the  attractor  well  enough  to  allow  it  to  be 
sufficiently  embedded  in  a  2d+l  phase  space. 

Certain  types  of  noise  make  it  possible  to  mistake  a  data  set  with  a  finite 
correlation  dimension  as  having  an  attractor.  Certain  sets  of  stochastic  data 

can  yield  finite  correlation  dimensions  yet  arc  not  strange  attractors.  For 
instance,  the  “random  walk”  yields  a  correlation  dimension  of  1.1  (Ramsey  and 
Yuan,  1989).  Even  more  disconcerting  (especially  to  those  meteorologists  who 
live  and  die  by  the  -5/3  power  law)  is  that  “colored”  random  noise 
characterized  by  a  power  law  spectrum  can  yield  a  finite  correlation 
dimension  (Osborne  and  Provenzale,  1989).  White  noise  which  has  a  flat  power 

spectrum  does  yield  an  infinite  correlation  dimension  as  indicated  originally 
by  Grassberger  and  Procaccia  (1983,  1984).  Osborne  and  Provenzale  (1989) 
showed  that  white  noise  gives  infinite  correlation  dimensions  because  the 
random  noise  acts  as  a  fractal  path  in  phase  space,  leading  to  self  similarity. 

3.2.2  Filtering  and  Digitizing  Errors 

Filtering  and  digitizing  the  data  can  yield  inaccurate  estimates  of  the 
dimension  of  an  attractor.  Filtering  leads  to  inaccurate  dimension  estimates  in 

much  the  same  manner  as  noise.  Data  is  often  filtered  to  reduce  the  amount  of 
noise  in  a  signal,  but  doing  so  can  lead  to  an  overestimate  of  the  dimension  of 
an  attractor  (Badii  ct  al.,  1988).  On  the  other  hand,  digitizing  the  signal  creates 
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errors  which  underestimate  the  dimension;  in  this  regard  it  acts  in  an  opposite 
manner  to  noise  (Mdller  et  al.,  1987).  Random  noise  combined  with  errors 
introduced  during  digitizing  can  create  a  "false"  scaling  region  which  could  in 
turn  suggest  the  presence  of  an  attractor  where  there  was  none. 

Mdller  et  al.,  (1987)  suggest  that  adding  Gaussian  noise  with  a  standard 
deviation  equal  to  0.4  times  the  least  significant  bit  before  digitizing  reduces 
the  error  in  the  dimension  estimate  in  signals  where  most  of  the  error  is 
attributed  to  digitizing.  For  signals  primarily  distorted  by  noise,  adding 
Gaussian  noise  with  a  standard  deviation  on  the  order  of  the  least  significant 
bit  provides  the  best  improvement  in  the  dimension  estimate.  Tests  run  by 
Mdllcr  et  al.,  (1987)  for  both  cases  show  the  error  in  the  dimension  estimate 
can  be  reduced  by  a  factor  up  to  80%. 

3.2.3  Number  of  Points 

A  limited  data  series  (i.e.  limited  in  the  total  number  of  points  or  amount 
of  the  attractor  that  is  covered)  leads  to  a  downward  bias  in  the  dimension  of 
random  variables  and  an  upward  bias  in  the  the  estimate  of  a  dimension  of  an 
attractor  (Ramsey  and  Yuan,  1989).  Small  data  seu  also  lead  to  conditions  where 
the  correlation  integral  does  not  saturate  at  increasing  embedding  dimension. 
Ramsey  and  Yuan  (1989)  suggest  a  method  of  non-linear  curve  fitting  that  will 
allow  one  to  test  for  the  presence  of  an  attractor  in  a  limited  data  set. 

Obviously,  there  must  be  some  minimum  number  of  points  for  which 
the  Grassberger-Procaccia  algorithm  will  yield  accurate  estimates  of  the 
attractor  dimension.  The  most  commonly  quoted  limits  on  the  number  of  points 
to  adequately  implement  this  algorithm  is 


10J2-  I00d2 


(3.26) 
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where  d2  is  the  correlation  dimension  (see,  for  example,  Henderson  and  Wells, 
1988). 

Smith  (1989)  provides  a  more  rigorous  and  detailed  formulation  of  the 
minimum  number  of  points  needed  to  obtain  an  accurate  dimension  estimate. 
Smith  gives  the  number  of  points  necessary  to  estimate  the  correlation 
dimension  of  a  non-lacunar  set  to  within  5%  of  its  true  dimension  as 


Nmh>42M 

(3.27) 

where  M  is  the  greatest  integer  less  than  the  dimension.  For  example,  the 
value  of  M  would  be  2  for  an  attractor  with  a  dimension  of  2.3. 

Abraham  ct  al.  (1988)  showed  promising  results  in  examining  the 
dimensions  of  small  data  sets  contaminated  by  noise.  They  were  able  to  obtain 
the  dimension  of  the  Hcnon  attractor  for  data  sets  with  as  few  as  500  points. 
Note  that  this  fits  Smith's  minimum  criteria;  the  dimension  of  the  Henon 
attractor  is  1.24  -  hence  the  minimum  number  of  points  needed  to  estimate  its 
dimension  could  be  as  low  as  42!  While  it  was  difficult  to  accurately  determine 
the  dimension  from  small  data  sets,  Abraham  et  al.  (1988)  still  felt  it  was 
possible  to  distinguish  between  chaotic,  periodic  and  random  behavior. 

3.2.4  Time..  Delay 

A  single  time  scries  may  not  properly  fill  out  phase  space  if  the  wrong 
time  delay  for  the  embedding  dimension  is  chosen.  Often,  the  time  to  the  first 
zero  in  the  autocorrelation  is  chosen  as  the  delay  time  in  constructing  the 
higher  dimension  vectors  from  the  time  series.  Fraser  and  Swinney  (1986) 
points  out  this  practice  as  being  “naive”;  it  may  grossly  underestimate  the 
correlation  dimension  of  the  attractor  in  the  data.  Too  small  a  time  delay  in  the 
case  of  highly  autocorrelated  data  may  yield  pairs  of  points  that  lie  close 
together  because  they  are  closely  related  in  time  rather  than  their  lying 
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"accidentally"  close  together  because  they  are  on  the  attractor.  Highly 
autocorrelatcd  data  restricts  the  trajectories  of  a  signal  from  "filling"  out 
phase  space,  thus  restricting  the  information  that  may  be  gained  by 
examining  the  distances  between  them.  Figure  3.3  shows  an  example  of  a 
function  with  two  choices  of  time  delay,  one  of  which  fills  the  phase  space  and 
one  of  which  does  not. 

Fraser  and  Swinney  (1986)  suggest  that  a  better  choice  of  the  proper 
time  delay  can  be  made  on  the  basis  of  mutual  information  theory.  The  mutual 
information,  I,  is  defined  as 


1  = 


P(  X,Y)  log2 


'  P( X, Y)  ' 
P(X)P(Y) 


dXd  Y 


(3.28) 


where  X  is  the  measurement  at  time  t,  Y  is  the  measurement  at  time  t+i,  P(X.Y) 
is  joint  probability  density  and  P(X)  and  P(Y)  arc  the  respective  X  and  Y 
probability  densities.  For  logarithms  taken  to  base  2,  the  units  of  the  mutual 
information  is  in  bits.  The  mutual  information  measures  the  relationship 
between  two  signals  in  a  more  general  manner  than  the  autocorrelation, 
which  measures  the  linear  dependence.  Normally,  the  first  zero  in  the 
autocorrelation  between  two  data  vectors  implies  that  the  two  are  linearly 
independent.  However,  data  characterized  by  a  strange  attractor  are  usually 
highly  nonlinear,  thus  making  the  first  zero  in  the  autocorrelation  a  poor 
choice  for  the  time  delay.  For  nonlinearly  related  data,  Fraser  and  Swinney 
suggest  that  first  local  minimum  in  the  mutual  information  provides  the  best 
choice  of  time  delay  in  construction  of  higher  dimension  data  vectors. 
However,  calculating  the  mutual  information  is  computationally  very 
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expensive  (Fraser  and  Swinney,  1986)  and  often  not  as  enlightening  as 
repeating  the  correlation  integral  calculations  for  different  time  delays4. 

Liebart  and  Schuster  (1989)  show  the  first  local  minimum  in  the  mutual 
information  is  the  best  choice  for  time  delay  as  opposed  to  some  other  local 
minimum  (say  the  second  or  third).  Their  tests  show  that  the  first  local 
minimum  in  the  mutual  information  helps  preserve  the  small  scale  structure 
of  the  attractor  in  the  phase  space  reconstruction.  They  also  point  out  that  this 
criterion  for  the  time  delay  is  not  "that  the  reconstructed  orbit  in  phase  space 
is  closest  to  the  true  one  but  that  the  dimensions  and  entropies  from  the 
reconstructed  orbit  are  closest  to  their  true  values''.^ 

3.2.5  The  Norm 

Additional  error  can  be  introduced  by  an  improper  choice  of  norm  used 
in  the  Grassbergcr-Procaccia  method.  While  all  norms  are  theoretically 
equivalent,  Havstad  and  Ehlers  (1989)  found  differences  between  the 
dimensions  calculated  from  using  the  Euclidean  norm  and  the  maximum  notm. 
The  Euclidean  norm  counts  points  that  fall  within  spheres  while  the  maximum 
norm  counts  points  that  fall  within  cubes.  Havstad  and  Ehlers  (1989)  found  the 
maximum  norm  underestimated  the  dimension  of  the  Mackey-Glass  attractor 
whereas  the  Euclidean  norm  yielded  a  value  for  the  dimension  that  was  very 
close  to  the  true  one.  They  believe  that  the  difference  occurs  because  the 
diagonals  of  the  cubes  are  aligned  with  the  surfaces  of  the  attractor.  While  the 
maximum  norm  is  attractive  because  if  js  computationally  less  intensive,  only 
the  more  computationally  expensive  Euclidean  norm  will  be  used  in  this  study. 


4  Glenn  James,  personal  communication  (1989).  Glenn  James  is  a  fellow  AF 
Ph.D.  type  who  just  recently  graduated  from  Georgia  Tech  (Winter  quarter, 
1990).  He  also  noted  that  the  mutual  information  calculations  were  too  much 
trouble  and  "a  pain  in  the 
^Liebart  and  Schuster  (1989),  pl08. 


53 


3.2.6  Non-Stationarv  Data  Sets 

Implicit  in  this  analysis  is  the  assumption  that  the  data  set  is  stationary, 
an  assumption  common  to  most  signal  analysis.  However,  Havstad  and  Ehlers 
(1989)  have  shown  the  Grassbcrger-Procaccia  method  is  acceptable  (with  some 
modifications)  for  data  sets  that  are  not  stationary.  The  dimension  for  small, 
non-stationary  data  sets  can  be  calculated  from  small  overlapping  groups  with 
reasonable  accuracy  if  the  dimension  is  not  more  than  10  (Havstad  and  Ehlers, 
1989).  However,  the  number  of  points  in  each  overlapping  segment  still  must 
be  the  minimum  number  necessary  to  implement  the  correlation  integral 
algorithm. 

Given  some  of  the  errors  and  uncertainties  discussed  above  (small  data 
sets,  noise,  a  certain  degree  of  autocorrelation,  non-stationary  data  sets),  one 
of  the  most  difficult  problems  is  defining  the  scaling  region.  In  many 
experimental  results,  the  scaling  region  is  quite  small  and  poorly  defined  for 
limited  sets  of  noisy  data.  Ellner  (1988)  has  developed  an  alternate  method  of 
calculating  the  generalized  dimension  of  an  attractor  that  improves  the 
discrimination  of  the  scaling  region  as  well  as  providing  an  estimate  of  the 
errors.  This  technique  is  based  on  a  maximum  likelihood  method. 

Unfortunately,  the  maximum  likelihood  method  is  even  more  computationally 
demanding  than  the  Grassberger-Procaccia  algorithm. 

Ellner  (1988)  notes  that  the  maximum  likelihood  technique  offers 
several  advantages  over  the  technique  of  Grassbcrger  and  Procaccia.  For  small 
data  sets,  the  scaling  region  is  better  defined  and  freer  of  distortions  that 
occur  because  of  the  finite  sample  effect.  This  yields  a  larger  apparent  scaling 
region  and,  in  turn,  gives  a  more  accurate  estimate  of  the  dimension. 
Furthermore,  the  maximum  likelihood  technique  yields  a  dimension  estimate 
accompanied  by  confidence  intervals  which  give  the  error  due  to  finite 
sample  size. 
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3.3  Extension  of  the  Grassbereer-Procaccia  Correlation  Integral 


Additional  research  has  been  done  examining  the  strengths  and 
weaknesses  of  this  technique.  One  of  the  biggest  disadvantages  of  the 
technique  proposed  by  Grassberger  and  Procaccia  (1983)  is  the  limits  in 
equations  3.15  and  3.18.  These  limits  impose  a  requirement  for  a  lengthy  time 
series  and,  as  pointed  out  by  Theiler  (1986),  can  lead  to  spurious  dimensions  if 
the  number  of  data  points  is  too  small  and  the  data  are  too  highly 
autocorrelated.  However,  Theiler  proposes  a  modification  to  the  technique  of 
Grassberger  and  Procaccia  (1983)  which  improves  the  convergence  of  the 
integral  towards  its  infinite  limit  for  autocorrelated  data. 

Theiler  redefines  the  correlation  integral  of  Grassberger  and  Procaccia 
(1983)  as: 


C(e,N)  =  -2- £  £  H(£-|xi+n-Xj|) 
n=l  i=l 

(3.29) 

where  H(x)  is  the  Heaviside  step  function,  N  the  total  number  of  points  in  the 
data  set  and  the  phase  space  vector  xj  is  defined  as 

Xj  =  (sj,  Sj+t,  Si+2x,  •  •  •,  Si^m.i)T)  where  sj  is  the  original  signal  at  time  t  and  t  the 

time  delay.  Theiler  notes  that  for  typical  conditions  (t,  m  «  N)  .there  are  as 
almost  as  many  vectors  xj  as  there  are  data  points  sj.  The  correlation 
dimension,  v,  can  be  defined  by  the  limit 


v=  lim  lim 

e— >0  N  — >oo  e 


(3.30) 


or  where  the  derivative  exists 
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v  =  lim  lim 
£— >0  N  — >oo 


d[  log  C(e,  N)] 

_ d£ _ 

dloge 

de 

(3.31) 

The  limits  in  the  two  expressions  above  provide  problems  in  implementing 
this  technique  in  practice. 

Theiler  (1986)  shows  that  the  Grassberger-Procaccia  correlation 
integral  redefined  in  equation  3.29  can  be  better  served  by  a  more  generalized 
version 


N  N-n 

C(e,N,  W)  =  -\  X  X  H(E-|xi+„-Xi|) 
N  n=Wi=l 


(3.32) 


W  is  the  number  of  autocorrelated  points  to  exclude  from  the  counting 
statistics.  Note  that  W=1  yields  the  original  definition  of  the  correlation 
integral  proposed  by  Grassbergcr  and  Procaccia  (1983).  Essentially,  this 
algorithm  doesn't  count  the  first  W  autocorrelated  points  that  lie  nearby. 
Skipping  these  points  in  the  summation  improves  the  convergence  of  the 
algorithm  in  both  limits.  Theiler  notes  that  the  key  to  this  modification  lies  in 
choosing  the  right  value  of  W. 

Theiler  (1986)  demonstrated  the  validity  of  this  modification  in 
examining  uncorrelated  and  autocorrelated  noise.  As  previously  mentioned, 
autocorrelation  restricts  the  trajectories  of  a  function  from  “filling”  out  phase 
space,  thus  restricting  the  information  that  may  be  gained  from  examining 
the  distances  between  pairs  of  data  points.  Highly  autocorrelated  data  may 
yield  pairs  of  points  that  lie  close  together  because  they  are  closely  related  in 
time  rather  than  their  lying  “accidently”  close  together  because  they  are  on 


an  attractor. 


Noise  has  a  slope  on  the  C(e,N)  vs  e  plots  equal  to  m  since  noise  scales 
as£m  where  m  is  the  embedding  dimension  (Theiler,  1986),  Noise  has  an  effect 
similar  to  that  for  autocorrelation.  Autocorrelation  and  a  limited  number  of 
points  may  unnecessarily  restrict  the  range  over  which  this  slope  occurs, 
creating  a  "knee"  in  the  plot  of  In  C(e,N)  vs.  In  e  at  higher  embedding 
dimensions.  Theiler  (1986)  shows  that  the  usable  range  of  C(e,N,W),  i.e.  that 
range  over  which  C(e,N,W)  is  proportional  to£m,  will  be  between 

a.  Uncorrelated  limit  -  2/N2  and  1,  and 

b.  Autocorrelated  limit  -  2/N2  and  2/N. 

As  the  number  of  data  points  approaches  infinity,  the  autocorrelated  range 
may  approach  the  uncorrelated  limit. 

Given  this,  Theiler  suggests  a  minimum  value  for  W  which  extends  the 
usable  range  for  the  slope: 


w  >t(^)2/m 

(3.33) 

where  x  is  the  first  zero  in  the  autocorrelation,  N  the  number  of  data  points 
and  m  is  the  embedding  dimension.  This  is  equivalent  to  dropping  the  first  W 
terms  in  the  summation  series  in  the  correlation  integral.  By  dropping  the 
first  W  terms,  the  summation  neglects  the  points  that  are  nearby  because  they 
are  correlated  in  time  and  thus  approaches  its  true  limit. 

There  are  two  different  time  scales  which  must  be  considered  in 
employing  the  Grassberger-Procaccia  algorithm.  The  time  scale  of  the  first 
minimum  of  the  mutual  information  determines  the  best  lag  for  the 
reconstruction  of  phase  space  vectors  from  the  original  data.  The  time  to  the 
first  zero  in  the  autocorrelation  determines  the  number  of  autocorrelated  data 
points  to  exclude  from  the  summation  in  the  correlation  integral. 
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3.3.1  The  Normalization  Factor 


There  is  some  question  about  the  normalization  factor  used  in  the 
correlation  integral  algorithm  and  its  variations.  In  the  original  definition  by 
Grassberger  and  Procaccia  (1983,  1984),  the  normalization  factor  was  1/N2. 
However,  it  was  seen  that  the  number  of  calculations  could  be  halved  by  only 
determining  the  distances  for  the  upper  half  of  the  matrix  which  held  the 

pairs  of  points  (i.e.,  the  distance  between  point  i  and  point  j  is  the  same  as  the 

distance  between  point  j  and  point  i  —  why  repeat  the  calculation?).  The 
normalization  factor  was  further  modified  by  Henderson  and  Wells  (1988)  to 
2/N(N-l)  by  not  calculating  the  identity  (sometimes  called  "self  pair")  terms 

(i.e.,  those  terms  for  which  i=j).  Smith  (1988)  noted  that  the  "self  pair"  terms 

must  be  calculated  explicitly  and  included  in  the  full  double  summation.  This  is 
required  to  distinguish  the  scaling  of  true  noise  from  fluctuations  due  to  a 
finite  number  of  points  N.  When  the  "self  pair"  (i=j)  terms  are  omitted,  the 
correlation  integral  is  not  necessarily  bounded  in  the  limit  as  e  approaches 
zero. 


3.4  Autocorrelation  Analysis 

The  autocorrelation  analysis  is  based  on  the  correlation  function  which 
is  defined  as  follows  (Walpole  and  Meyers,  1989): 


r  = 


1  x  y 


V  SXX  Syy 


where  Sxx,  Syy  and  SXy  are  given  by 


(3.34) 
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The  individual  components  in  the  summation  are  defined  as 

x;  =  x(iAt) 
yi  =  x((i  +  X)At) 

(3.36) 

where  xj  is  an  element  of  the  original  time  series  x.  At  is  the  sampling  interval 
and  yj  is  an  element  of  the  original  time  series  x  lagged  by  some  factor  X. 

The  correlation  was  calculated  for  each  lag  X.  The  error  of  the 
autocorrelation  is  defined  as 


)_ 

N-X 

1  - 

N(N  +  2) 

(3.37) 


The  errors  were  calculated  but  are  not  displayed  on  the  figures  shown  in  this 
study  simply  to  render  the  figures  more  readable.  The  number  of  points  in  the 
autocorrelation  were  chosen  to  minimize  the  error. 
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3.5  Power  Spectrum  Analysis 


We  will  use  the  method  of  calculating  the  power  spectrum  as  described 
by  Press  et  al„  (1986).  The  "power"  calculated  in  this  study  is  the  mean  squared 
power  defined  as 


(3.38) 


where  A  is  the  sampling  interval,  N  the  number  of  points  and  T=(N-1)A.  The 
power  (P)  is  defined  as  a  function  of  frequency  (f)  at  N/2+1  discrete 
frequencies  by 


P(0)  =  -L-|coP 

p|  fk) = -y 1 ck|2 + 1  cN-k,2J ;  k=i' 2 . (i  -  >) 

p(fcl  =  -^|cwJ2 

where  the  frequencies  f^  and  fc  are  defined  as 


(3.39) 


f  =  — 

1  C 

2A 

fk  =  _k_  =  2fcJt  ;  k=0,  1, 

NA  N 


N 

2 


and  the  Fourier  coefficients  cK  arc  defined  by  the  following 


(3.40) 
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N-l 

ck  =  Icje^jk/N;  k=0,  1 . N-l 

j=o 


(3.41) 


The  "Welch"  window  was  used  on  the  data  when  applying  the  Fourier 
transform  (Press  et  al.,  1986).  The  data  was  segmented  to  obtain  the  smallest 
variance  of  the  power  for  the  number  of  data  points  in  the  transform.  The 
segments  were  overlapped  by  half  their  length  M  and  the  variance  was 
reduced  by  a  factor  of  9k/ll  where  k  is  the  number  of  segments.  The  number 
of  data  points  required  by  the  transform  is  (2k+l)M. 


61 


Figure  3.1  Correlation  integral  of  the  sine  function  with  period  equal  to 
thirty  with  10%  external  noise.  The  total  number  of  points  is 


SEE 


Kolmogorov  entropy  (K2)  for  the  correlation  integral  of  the  sine 
function  shown  in  Figure  3.1  plotted  as  a  function  of  increasing 
embedding  dimension.  The  curve  asymptotically  approaches  the 
value  of  0,  indicating  periodic  behavior. 
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Figure  3.3  Phase  space  portrait  of  a  sine  function  with  period  of  300  for  a 
time  delay  of  x~\  and  a  time  delay  equal  to  the  first  zero  in  the 
autocorrelation,  i.e.  x=75. 


CHAPTER  IV 


OVERVIEW  OF  THE  DATA 


This  chapter  gives  an  overview  of  the  data  used  in  this  study.  It  first 
describes  the  theory  behind  the  partial  reflection  radar  from  which  the  data 
was  obtained.  This  is  followed  by  a  description  of  the  Saskatoon  partial 
reflection  radar  itself.  The  chapter  concludes  with  a  description  of  the  data 
from  the  Saskatoon  radar  used  in  this  study. 


4.1  Theory  of  Partial  Reflection  Radars 

Partial  reflection  radars  transmit  an  electro-magnetic  pulse  which  is 
both  partially  and  totally  reflected  in  the  middle  atmosphere.  The  reflected 
signal  creates  a  diffraction  pattern  on  the  ground  which  moves  at  a  rate  twice 
the  speed  of  the  wind  in  the  scattering  layer  at  which  the  reflection  occurred 
(see  Figure  4.1).  This  induces  an  e.m.f.  in  a  stationary  antenna.  As  the  pattern 
moves  past  the  antenna,  the  induced  e.m.f.  varies  in  amplitude;  this  is  called 
radiowave  fading.  Normal  signal  analysis  of  partial  reflection  radar  data 
revolves  around  cross  correlation  of  the  signal  peaks  between  antennas  to 
deduce  the  direction  and  speed  of  the  wind.  Comparisons  of  the  fading  rates 
from  at  least  three  antennas  give  an  estimate  of  the  horizontal  pattern 
velocity  and  hence  the  velocity  of  the  reflecting/scattering  region. 

The  mechanisms  causing  the  echoes  range  from  reflections  from 
sharply  bounded  irregularities  in  the  refractive  index,  to  scattering  (i.e., 
"partial  reflections")  from  quasi-isotropic  irregularities  ("turbulent  blobs") 
whose  scales  are  on  the  order  of  half  the  radar  wavelength.  At  the  altitudes 
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measured  by  partial  reflection  radars  (typically  50  -  120  km,  depending  on  the 
design  and  sensitivity  of  the  radar),  irregularities  in  the  refractive  index  can 
be  created  by  free  electrons  carried  along  by  the  neutral  wind  as  well  as 
variations  in  potential  temperature  and  moisture  content  (Hocking,  1985). 

Short  wavelength  solar  radiation  ionizes  the  air  in  the  upper  atmosphere  and 
generates  free  electrons;  concentrations  of  these  free  electrons  typically  vary 
from  less  than  ~  102  cm-3  at  60  km  to  ~  104  cm“3  at  100  km  (Gregory  and 
Stephenson,  1972). 

The  refractive  index  varies  primarily  in  the  vertical,  a  consequence  of 
the  increase  in  the  free  electron  density  with  height  as  well  as  changes  in  the 

potential  temperature  with  height.  The  signal  strength  generally  increases 
with  height  because  of  the  vertical  gradient  in  the  refractive  index.  However, 

there  are  horizontal  variations  in  the  refractive  index  which  modulate  the 
transmitted  pulse  wave  front  so  that  the  scattered/reflected  wave  is  equivalent 

to  an  angular  spectrum  of  plane  waves.  The  horizontal  variations  in  the 

refractive  index  can  be  either  a  result  of  horizontal  variations  in  the  potential 
temperature  caused  by  turbulence  or  wave  activity  or  horizontal  variations  in 
the  free  electron  density  or  a  combination  of  both.  This  spectrum  of  plane 
waves  creates  the  diffraction  pattern  on  the  ground  which  is  measured  by  the 
receiving  antennas  (Fraser,  1984), 

Echoes  can  be  obtained  from  minimum  heights  of  50-60  km  to  heights 
where  the  signal  is  totally  reflected  in  the  E  or  F  region  of  the  ionosphere1. 

The  lower  altitude  limit  is  essentially  a  function  of  the  sensitivity  of  the  system 
to  signals  that  are  only  partially  reflected  in  the  lower  mesosphere.  The  spatial 
characteristics  of  the  partial  reflection  regions  vary  from  thin,  stratified 
layer_  to  thick,  turbulent  layers.  The  thin,  stratified  layers  are  often  less  than 
1  km  thick  (Fraser,  1984).  Measurements  made  near  local  solar  noon  can  be 
influenced  by  the  D  region  of  the  ionosphere  which  extends  further  down  into 

*The  E  region  spans  the  altitude  range  90  -  150  km.  The  F  region  encompasses 
150  -  500  km  (Kelley,  1989). 
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the  mesosphere2.  The  effects  of  the  geomagnetic  field  on  the  motions  of  the 
free  electrons  must  be  considered  at  heights  where  the  electron-neutral 
collision  frequency  is  low.  This  effect  does  not  appear  to  be  significant  at 
heights  below  95-100  km  (Fraser,  1984). 


4.2  Saskatoon  Partial  Reflection  Radar 

The  data  used  in  this  study  came  from  the  partial  reflection  radar 
located  at  52°  N,  107°  W  in  Saskatoon,  Canada.  This  radar  operates  at  an  average 
frequency  of  2.2  MHz  with  an  equivalent  wavelength  of  135  m.  It  emits  an 
approximately  trapezoidal  pulse  with  a  width  of  20  ps;  this  is  equivalent  to  a  3 
km  height  resolution  (Manson  et  al.,  1974;  Gregory  and  Stephenson,  19'72; 
Manson  and  Meek,  1987;  Meek  ard  Manson,  1987).  Data  are  recorded  for  23 
height  levels  ranging  from  52  km  to  118  km  in  3  km  increments.  The  pulse 
repetition  rate  used  for  this  data  set  is  15  s-1.  Data  below  70  km  is  considered 
unreliable  (Meek.  1989,  private  communication). 

The  radar  consists  of  a  transmitter  and  four  receiving  antennas.  The 
receiving  antenna  array  is  laid  out  in  a  “Y”  pattern  with  a  separation  between 
antennas  1,  2  and  3  of  2\  (270  m),  twice  the  wavelength  of  the  radar.  The 
receiving  antenna  array  is  depicted  in  Figure  4.2.  Measurements  are  taken  at 
each  of  the  four  antennas  as  the  recording  system  cycles  around  the  antenna 
array  at  a  rate  of  15  Hz;  this  yields  the  0.2666  second  measurement  separation 
between  measurements  at  each  antenna. 

4.3  Overview  of  the  Data 

The  data  set  was  taken  on  2  August  (day  of  year  214)  1985  starting  at 
18:31:00  GMT  and  covers  approximately  one  and  a  half  hours  (a  total  of  17910 

2Thc  D  region  is  less  than  90  km  altitude. 
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points).  The  data  set  has  been  described  previously  in  the  literature  by  Manson 
and  Meek  (1987).  Manson  and  Meek  (1987)  found  several  "glints"  or  regions  of 
focused  turbulence  in  the  data  set  to  be  used  in  this  study. 

The  data  were  brokr  down  into  three  separate  files.  There  is  a  time  gap 
between  the  files  due  to  data  recording  considerations.  The  first  file  started  at 
18:31:00  GMT,  ended  at  19:01:12  GMT  and  consisted  of  6705  points.  There  was  a  2 
minute  16  second  gap  between  the  first  file  and  the  second  file  which  is 
equivalent  to  510  points.  The  second  file  began  at  19:03:28  GMT  and  ended  at 
19:29:36  GMT  for  a  total  of  5895  points.  There  is  a  gap  of  approximately  40 
seconds  between  the  second  and  the  third  files  which  is  equivalent  to  146 
points.  The  third  file  began  at  19:30:15  GMT  and  ended  at  19:53:50  GMT  for  a  total 
of  5310  points.  The  total  number  of  points,  not  counting  the  breaks  between 
files,  is  17910. 

A  careful  examination  of  the  raw  data  record  seemed  to  reveal  an 
undocumented  change  of  gain  between  the  first  and  remaining  files.  The 
apparent  change  in  gain  was  signalled  by  a  20  second  dropout  in  signal  just 
prior  to  the  end  of  the  first  file  (i.e.,  at  point  6610).  We  calculated  the  variance 
of  the  data  for  the  first  6000  points  and  com;  ared  that  to  the  variance  of  the 
remaining  points  for  the  seven  levels  between  76  and  94  km.  The  results  of 
these  calculations  are  summarized  in  Table  4.1.  The  average  ratio  of  the 
variances  over  all  seven  levels  was  3.08  which  is  very  close  to  what  you  would 
expect  from  a  10  db  (i.e.,  a  ratio  of  3.16)  change  in  gain. 

We  used  the  F  test  (Walpole  and  Meyers,  1989)  to  determine  whether  the 
variances  were  indeed  the  same.  We  were  able  to  reject  the  hypothesis  that  the 
variances  were  the  same  for  each  level  at  the  99.995%  confidence  level.  Thus, 
it  is  safe  to  conclude  that  there  was  an  undocumented  gain  change  between 
the  first  and  the  remaining  files. 
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Table  4.1  Summary  of  mean  and  variance  changes  in  the  data. 


Level  (km) 

1  <  n  <  6000 

Mean  Variance 

7000  <  n  <  17900 

Mean  Variance 

Ratio  of  two 

variances 

76 

101.4 

306.4 

103.8 

99.8 

3.07 

79 

100.4 

969.7 

104.4 

680.3 

1.43 

82 

99.0 

1569.8 

102.6 

401.3 

3.91 

85 

100.3 

1052.6 

102.9 

413.3 

2.55 

88 

98.1 

1024.2 

102.3 

288.7 

3.55 

91 

111.8 

2341.9 

108.4 

400.9 

5.84 

94 

1 15.4 

958.0 

113.9 

774.8 

1.24 

There  is  also  a  difference  in  the  mean  signal  strength  for  each  of  the 
seven  levels  between  the  first  file  and  the  last  two  files.  We  used  the  Student  t 
test  to  determine  whether  the  means  belonged  to  the  same  population.  We  were 
able  to  reject  the  hypothesis  that  the  mean  from  the  first  file  and  the  mean 
from  the  second  and  third  files  at  each  level  belonged  to  the  same  population 
at  the  99.95%  level.  The  change  in  mean  over  the  approximate  hour  and  a  half 
that  is  spanned  by  the  data  set  is  most  likely  due  to  receiver  drift.  Meek 
(personal  communication,  1989)  reports  that  the  data  should  be  corrected  for 
receiver  drift  for  observation  periods  longer  than  an  hour. 

The  variance  is  low  for  the  76  km  data  compared  to  that  for  the  other 

levels  above  76  km.  In  general,  the  signal  strength  increases  with  height. 

There  was  very  little  variation  in  the  signal  for  levels  below  76  km  which  is 

why  data  from  these  levels  were  not  even  considered  beyond  some 

preliminary  examination. 

We  tested  for  receiver  drift  over  the  lengtn  of  the  first  file  by 
calculating  the  mean  and  variance  in  500  point  groups  from  point  number  one 
to  point  number  6500  for  each  of  the  seven  levels  (76  through  94  km).  The 
hypothesis  that  the  means  for  each  level  belonged  to  the  same  population  was 
accepted  at  the  95%  confidence  level.  While  this  hypothesis  was  not  accepted  at 
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a  more  rigorous  confidence  level  (i.e.,  one  >  95%),  it  seems  reasonable  to 
assume  that  any  receiver  drift  over  the  first  data  file  (a  period  of 
approximately  30  minutes)  was  minimal. 

These  two  problems,  receiver  drift  and  change  in  signal  gain,  could  be 
corrected..  The  drift  in  the  mean  could  be  easily  remedied  by  detrending  the 
data.  The  change  in  the  gain  could  be  eliminated  by  subtracting  out  the 
detrended  mean,  multiplying  by  the  ratio  of  the  variances  and  adding  back  the 
mean.  While  making  any  changes  to  the  original  data  increases  the  risk  of 
rendering  any  dimension  estimates  inaccurate,  these  changes  represent  the 
least  significant  risk.  However,  there  was  one  additional  problem  that 
eliminated  almost  two  thirds  of  the  data  from  consideration. 

In  addition  to  the  gain  change,  the  data  in  the  second  and  third  files 
appeared  to  be  contaminated  by  a  high  frequency  oscillation.  This  was  almost 
certainly  due  to  a  folding  back  of  some  portion  of  the  long  range  signal  into 
the  lower  height  gates.  The  oscillation  can  be  seen  clearly  in  the 
autocorrelation  (see  Figure  4.3)  and  power  spectrum  (see  Figure  4.4)  of  data 
from  the  second  and  third  files.  The  oscillation  was  present  in  the  signal  at  all 
of  the  heights  (76-94  km)  to  be  used  in  this  study,  but  it  was  strongest  in  the 
data  from  91  and  94  km.  The  oscillation  occurred  with  a  period  of  .867  s,  which 
is  slightly  greater  than  every  third  point. 

Filtering  the  data  could  easily  remove  this  high  frequency  signal. 
However,  doing  so  would  increase  the  uncertainty  in  any  estimate  of  the 
dimension  of  the  system  and  possibly  suggest  the  presence  of  an  attractor 
where  there  was  none.  Because  of  this,  we  decided  to  eliminate  data  from  the 
second  and  third  data  files  from  further  analysis.  Only  data  from  the  first  file 
(6705  points  -  18:31:00  to  19:01:12  GMT)  was  considered  in  further  analysis. 

The  data  for  each  of  the  four  antennas  and  seven  levels  from  76  to  94  km 
is  shown  in  Figures  4.5  -  4.8.  The  data  shown  in  the  figures  was  averaged  over 
3.5  s  and  only  every  thirteenth  point  was  plotted  in  order  to  improve  legibility. 
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Several  features  are  immediately  obvious  from  examining  the  data.  The 
variation  in  the  signal  for  levels  76  and  79  km  is  generally  much  less  than  for 
the  upper  levels  for  all  antennas.  In  addition,  the  variation  in  signal  strength 
is  slightly  less  at  all  levels  for  antenna  #4  than  for  the  other  three  antennas. 

Meek  and  Manson  (1987)  found  a  strong  scattering  layer  to  be  located  at 
82  km  in  a  previous  study  of  this  data  set.  The  variation  in  signal  strength  then 
increases  again  for  91  and  94  km.  This  phenomena  can  be  most  easily  seen  in 
the  data  from  antenna  #4  shown  in  Figure  4.8.  The  presence  of  a  strong 
scattering  layer  near  82  km  is  often  seen  in  the  summer  mesosphere  (Fraser, 
1984). 
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Figure  4.2.  Receiving  array  of  Saskatoon  radar  in  the  Y  configuration, 
(adapted  from  Meek  and  Manson,  1987) 
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Figure  4.3  Autocorrelation  for  antenna  #1  at  94  km  for  points  14000 
through  16000  (i.e.,  data  from  the  contaminated  third  file). 


Figure  4.5  Data  from  antenna  #1  for  (a)  76  km;  (b)  79  km;  (c)  82  km;  (d)  85 
km;  (e)  88  km;  (0  91  km;  and  (g)  94  km.  The  data  starts  at  18:31:00 
GMT.  The  figure  is  on  the  following  pages. 
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Figure  4.5  (g) 


t 


Figure  4.6  Data  from  antenna  #2  for  (a)  76  km;  (b)  79  km;  (c)  82  km;  (d)  85 
km;  (e)  88  km;  (f)  91  km;  and  (g)  94  km.  The  data  starts  at  18:31:00 
GMT.  The  figure  is  on  the  following  pages. 
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Figure  4.6  (b) 
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Figure  4.6  (c) 
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Figure  4.6  (d) 
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Figure  4.7  Data  from  antenna  #3  for  (a)  76  km;  (b)  79  km;  (c)  82  km;  (d)  85 
km;  (e)  88  km;  (f)  91  km;  and  (g)  94  km.  The  data  starts  at  18:31:00 
GMT.  The  figure  is  on  the  following  pages. 
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Figure  4.8  Data  from  antenna  #4  for  (a)  76  km;  (b)  79  km;  (c)  82  km;  (d)  85 
km;  (e)  88  km;  (0  91  km;  and  (g)  94  km.  The  data  starts  at  18:31:00 
GMT..  The  figure  is  on  the  following  pages. 
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CHAPTER  V 


ANALYSIS  AND  INTERPRETATION 

This  chapter  contains  the  results  of  the  analysis  of  the  data  and  the 
interpretation  of  those  results.  It  first  describes  how  the  Grassberger- 
Procaccia  algorithm  was  applied  to  the  data,  followed  by  the  results  of  that 
analysis.  Analysis  of  the  power  spectra  of  the  echo  data  from  the  partial 
reflection  radar  follows.  The  chapter  concludes  with  the  interpretation  of 
results  from  the  implementation  of  the  Grassberger-Procaccia  algorithm. 


5.1  Analysis 

Determining  the  optimum  time  delay  for  building  the  phase  space 
reconstruction  of  the  attractor  is  the  first  step  in  applying  the  Grassberger- 
Procaccia  algorithm.  We  will  use  the  first  local  minimum  in  the  mutual 
information  as  the  best  time  delay  and  the  first  zero  in  the  autocorrelation  to 
determine  the  number  of  autocorrelated  vectors  to  exclude  from  the 
summation  in  the  correlation  integral,  i.e.,  the  value  of  W  from  equation  3.32. 

5.1.1  Calculating  the.  Jbest  choice  for  time  delay 

The  mutual  information  was  calculated  from  equation  3.28  with  up  to  200 
lags.  The  radar  data  ranges  from  a  minimum  value  of  0  to  a  maximum  value  of 
255;  the  integral  was  divided  into  boxes  with  sides  AX  and  AY  of  length  10.  This 
box  size  was  chosen  to  provide  the  best  convergence  for  the  integral  over  the 
entire  range  of  lags!  The  mutual  information  was  calculated  using  the  first 
6000  points  of  the  data  from  each  of  the  four  antennas.  It  was  not  possible  to 


calculate  the  mutual  information  for  a  smaller  number  of  points  because  the 
probabilities  would  not  be  statistically  valid.  An  even  larger  number  of  points 
would  have  been  desirable. 

The  results  of  the  mutual  information  calculation  are  depicted  in 
Figures  5.1  -  5.4.  In  general,  the  mutual  information  drops  off  very  quickly 
with  increasing  values  of  lag.  There  are  some  oscillatory  characteristics  at 
larger  lags  but  they  are  small  in  amplitude.  The  lag  for  which  the  first  minima 
in  the  mutual  information  occurs  is  summarized  for  each  level  in  Table  5.1;  the 
subscripts  indicate  the  antenna  for  which  the  calculations  were  performed. 


Table  5.1  Lac  for  the  first  local  minimum  in  the  mutual  information. 


Level  (km) 

*1 

n 

-C3 

U 

76 

12 

10 

9 

12 

79 

18 

20 

19 

18 

82 

17 

19 

18 

19 

85 

15 

19 

18 

15 

88 

14 

17 

16 

16 

91 

22 

19 

17 

17 

94 

13 

17 

14 

18 

5.1.2  A.UiQ.C.Orrg.l.ati.on 

The  autocorrelation  was  calculated  up  to  a  maximum  of  200  lags  using 
equation  3.35.  The  maximum  error  according  to  equation  3.37  was  less  than 
0.05%  at  the  two  hundredth  lag.  The  data  for  each  level  and  antenna  was 
divided  into  three  2000  point  groups  (1  <  n  <  2000,  2000  <  n  <  4000,  4000  < 

n  <  6000)  and  the  autocorrelation  was  calculated  separately  for  each  group. 

The  autocorrelations  for  each  antenna  and  level  are  shown  in  Figures 
5.5  -  5.8.  As  with  the  mutual  information,  the  autocorrelation  exhibited  some 


109 


oscillatory  behavior  as  a  function  of  lag.  In  general,  the  autocorrelation  was 
smoother  and  not  as  noisy  as  the  mutual  information. 

The  average  lag  (X)  at  which  the  first  zero  in  the  autocorrelation  occurs 
for  each  level  and  antenna  is  summarized  in  Table  5.2.  The  average  lag  is  the 
average  of  the  first  zero  of  the  three  2000  point  groups.  The  subscript  of  X 
indicates  the  antenna. 


Table  5.2  Lag 

(X)  for  the  : 

"irst  zero  in  t 

te  autocorrelation. 

Level  (km) 

Xi 

X2 

X3 

X4 

76 

24 

9 

8 

8 

79 

20 

18 

17 

20 

82 

22 

22 

28 

28 

85 

11 

10 

10 

11 

88 

15 

16 

13 

18 

91 

18 

18 

18 

18 

94 

15 

15 

17 

18 

5.1.3 


Since  we  intend  to  use  Theiler’s  modification  of  the  Grassberger- 
Procaccia  algorithm,  we  must  use  the  information  from  Tables  5.1  and  5.2  to 
calculate  the  value  of  W  (see  equation  3.33),  the  number  of  autocorrelated 
vectors  to  skip  in  the  summation.  The  value  of  W  for  each  antenna  and  level  is 
summarized  in  Table  5.3.  The  subscript  again  indicates  the  antenna.  The 
calculation  is  based  on  a  total  of  1500  points  and  a  maximum  embedding 
dimension  of  13.  The  values  have  been  rounded  up  to  the  next  greatest  integer 
since  Theiler's  modification  of  the  Grassberger-Procaccia  algorithm  employs 
integer  values  of  the  time  step  in  the  summation. 
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Table  5.3  The  value  of  W 

for  each  antenna  and  level 

N=1500,  m=13 

Level  (km) 

Wi 

w2 

w3 

W4 

76 

9 

4 

3 

3 

79 

8 

7 

7 

8 

82 

8 

8 

11 

11 

85 

4 

4 

4 

4 

88 

6 

6 

5 

7 

91 

7 

7 

7 

7 

94 

6 

6 

7 

7 

5.1.4  Implementing  the  Grassbereer-Procaccia  algorithm 

We  will  use  Theilcr's  (1986)  modification  of  the  Grassberger-Procaccia 
algorithm  to  calculate  the  base  statistics  necessary  to  estimate  the  dimension  of 
the  system.  The  exact  algorithm  implemented  in  this  analysis  is  given  by 

Jw" 

(5.1) 

where  N  is  the  number  of  points  in  the  data  series,  x  is  the  reconstructed  phase 
space  vector,  e  is  the  "search"  radius,  H  is  the  Heaviside  function  and  W  the 
number  of  autocorrelated  vectors  to  skip  (see  preceding  section  for  the  values 
of  W). 

We  used  the  normalization  factor,  (2/N(N-l)),  which  did  not  include  the 
{  calculation  of  the  "self  pair"  terms  nor  did  we  perform  the  summation  over  the 

full  range  of  indices.  Both  of  these  decisions  were  based  on  the  desire  to  reduce 
the  number  of  calculations  since  the  amount  of  computer  time  required  by  the 

algorithm  grows  by  the  square  of  the  number  of  points.  Tests  run  on  known 

r 

attractors,  strange  and  periodic,  and  noise  (i.e.,  the  Henon  attractor,  a  sine 
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function  of  known  period  and  Gaussian  noise)  with  a  limited  number  of  points 
(<1000)  revealed  no  convergence  problems  in  the  correlation  algorithm  as 
described  by  equation  5.1.  Since  the  Grassberger-Procaccia  algorithm  is  widely 
used  by  experimenters  without  calculating  the  full  double  summation  or  the 
"self  pair"  terms  and  our  tests  showed  no  decrease  in  convergence,  we  decided 
to  use  the  computationally  less  demanding  form  given  by  equation  5.1. 

The  algorithm  was  optimized  for  the  data  set.  The  correlation  integral 
was  calculated  for  150  values  of  e  evenly  spaced  logarithmically  over  the 
interval  from  In  e=l  to  In  e= 7.  This  was  repeated  for  12  different  values  of 
embedding  dimension  from  2  through  13.  Because  of  the  likelihood  of  the  data 
set  being  non-stationary,  the  algorithm  was  applied  to  overlapping  1500  point 
groups.  Each  group  overlapped  the  previous  group  by  half  the  number  of  data 
points  in  the  analysis.  Thus,  the  calculations  were  repeated  for  the  following 
seven  groups  of  points:.  1  <  n  <  1500;  750  <  n  <  2250;  1500  <  n  <  3000; 

2250  <  n  <  3750;  3000  <ni  4500;  3750  <  n  S  5250  and  4500  <  n  5  6000.  A 
limited  number  of  runs  were  done  for  5000  points.  Each  set  of  calculations  was 
repeated  for  each  level  (76  through  94  km)  of  each  antenna.  This  resulted  in 
over  196  separate  sets  of  calculations. 

All  the  calculations  were  performed  on  either  a  Macintosh  Ilex  or 
Macintosh  Ilci  computer.  FORTRAN  was  used  for  all  computer  programs. 

A  representative  set  of  plots  of  the  correlation  integral  for  each  level  is 
shown  in  Figure  5.9.  The  set  of  points  was  chosen  arbitrarily  as  was  the 
antenna.  This  set  of  figures  is  entirely  representative  of  all  the  overlapping 
segments  for  all  levels  for  each  of  the  antennas.  The  entire  set  of  calculations 
is  represented  by  the  graphs  shown  in  Appendix  A. 

5.2  Power  Spectrum  Analysis 

The  power  spectrum  was  implemented  upon  the  basis  of  equations  3.39- 
3.41.  The  Fast  Fourier  Transform  (FFT)  from  Press  et  al.  (1986)  required  the 


112 


number  of  points  analyzed  to  be  a  power  of  2.  To  obtain  a  power  spectrum 
which  extended  out  to  the  buoyancy  period  (~5  min)  required  a  total  of  8192 
points.  This  exceeded  the  total  number  of  points  in  the  first  file.  Consequently, 
some  data  (a  total  of  1487  points  --  approximately  6  min  36  s)  from  the  second 
file  were  included  in  the  power  spectrum  calculations.  We  did  not  correct  for 
the  gap  between  the  first  and  second  file. 

We  performed  tests  on  data  with  a  known  period  to  determine  the  degree 
of  aliasing  introduced  by  the  time  gap  between  the  files.  Tests  showed  that 
there  was  some  aliasing  but  also  that  it  was  not  expected  to  be  a  major  source  of 
error.  Any  errors  introduced  by  not  correcting  for  the  time  gap  are  on  the 
order  of  errors  introduced  by  the  change  in  variance  of  the  data  from  the 
second  file. 

The  power  spectra  for  antennas  1  through  4  at  all  levels  are  shown  in 
Figure  5.10-5.13,  respectively.  Note  that  the  x  axis  represents  the  period  and  is 
depicted  in  reverse  order.  This  was  done  to  give  the  figures  the  more  common 
appearance  of  power  spectra  as  a  function  of  frequency.  Interpreting  the 
power  spectra  as  a  function  of  period  is  the  same  for  frequency  as  long  as  one 
bears  in  mind  that  the  slope  of  any  power  law  behavior  will  be  the  negative  of 
what  you  would  normally  expect. 

The  power  spectra  exhibited  the  same  general  characteristics  for  all 
levels  and  antennas.  Each  power  spectrum  had  three  (sometimes  four)  distinct 
regions.  The  power  spectra  were  flat  for  periods  less  than  1-2  seconds.  In 
general,  the  spectra  were  flat  for  periods  less  than  2  seconds  for  76  km;  this 
region  decreased  to  periods  less  than  approximately  1  s  at  94  km.  The  flat 
region  was  followed  by  a  region  that  clearly  exhibited  power  law  behavior 
(i.e.,  constant  slope  on  a  log-log  plot);  this  usually  occurred  in  the  period 
range  between  1-2  and  approximately  6-10  s.  Weaker  power  law  behavior  (i.e., 
the  slope  on  the  log-log  plot  was  not  well  defined)  occurred  in  the  region 
between  approximately  10  and  200  s.  Each  power  spectrum  contained  the  peak 
energy  in  the  longest  period,  near  the  Brunt-Vaisdlla  period  for  the 


mesosphere.  Power  spectra  for  antenna  #4  were  qualitatively  similar  to  those 
for  the  other  three  antennas  except  that  the  power  over  the  entire  spectrum 
was  less.  This  is  most  likely  due  to  the  reduced  variance  in  the  data  from 
antenna  #4  described  in  Chapter  IV. 

To  further  understand  the  behavior  described  by  the  power  spectra,  we 
fit  different  portions  of  the  power  spectra  to  a  power  law.  These  results  are 
shown  in  Figures  5.14-5.17.  The  fitted  power  slopes  for  the  stronger  power  law 
behavior  (i.e.,  between  1-2  s  and  6-10  s)  are  summarized  in  Table  5.4. 


Table  5.4  Slo 

pe  for  power 

aw  curve  fits 

between  1-2  and  6-10  s. 

Level 

Antenna 

Antenna 

Antenna 

Antenna 

Average 

(km) 

#  1 

#2 

#3 

#4 

Slope 

76 

7.0 

5.3 

5.7 

5.2 

5.8 

79 

5.0 

6.6 

5.5 

- 

5.7 

82 

2.7 

2.3 

2.7 

1.9 

2.4 

85 

3.0 

2.9 

3.0 

2.4 

2.8 

88 

3.4 

1.5* 

3.1 

3.0 

2.8* 

91 

2.9 

1.5* 

3.5 

2.4 

2.6* 

94 

3.7 

2.1* 

4.2 

4.0 

3.5* 

*These  slopes  arc  distinctly  different  than  those  for  the  other  antennas.  The 
average  includes  these  values. 


We  must  be  careful  in  interpreting  the  apparent  power  law  behavior. 
There  has  been  almost  no  work  done  on  interpretation  of  power  spectrum 
analysis  of  echo  data  from  partial  reflection  radars  and  we  should  be  cautious 
in  treating  it  in  the  same  manner  as  that  for  horizontal  and  vertical  winds. 

Hocking  (1985)  defines  the  length  scale  for  the  viscous  region  in  the 
mesosphere  as  less  than  3-6  m  at  70  km  and  less  than  20-30  m  at  90  km.  In  the 
viscous  region,  the  kinetic  energy  of  turbulent  eddies  is  diminished  by  viscous 
effects  and  dissipated  as  heat. 
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The  length  scale  for  the  buoyancy  region  in  the  mesosphere  is  greater 
than  60-1000  m  (Hocking,  1985).  At  large  scales,  the  buoyancy  effects  become 
important.  In  the  buoyancy  region,  the  turbulent  eddies  take  on  a  "pancake" 
like  appearance  (Hocking,  1985)  and  consequently,  the  horizontal  length 
scales  are  much  greater  than  the  vertical  length  scales. 

The  inertial  region  lies  between  the  viscous  and  buoyancy  regions. 
Hocking  (1985)  notes  that  the  break  in  the  slope  of  power  spectra  between  the 
viscous  and  inertial  regions  occurs  at  length  scales  two  to  four  times  the 
length  scale  of  the  viscous  region. 

The  flat  region  in  the  power  spectra  corresponds  to  the  amount  of  time 
it  takes  the  smallest  size  eddy  that  can  be  detected  by  the  radar  (i.e., 
approximately  70  m)  to  be  advected  past  the  radar.  The  maximum  period  of  the 
flat  region  decreases  at  higher  levels  because  of  the  general  increase  in  the 
wind  speed  with  height.  Meek  and  Manson  (1987)  give  the  wind  speed  at  76  km 
as  approximately  30  m  s'1,  increasing  to  approximately  60  m  s*1  near  94  km 
for  this  time  period  at  Saskatoon. 

The  region  between  1-2  and  6-10  s  is  most  likely  the  viscous  region. 
While  the  radar  can  not  detect  turbulent  eddies  small  enough  to  be  in  the 
viscous  region,  it  can  detect  variations  in  the  motion  of  the  echo  patterns  that 
occur  on  time  scales  which  correspond  to  the  viscous  region.  The  power 
spectra  length  scales  for  the  viscous  region  can  be  calculated  from  the  period 
using  the  multiplicative  factor  from  Hocking  (1985),  the  wind  speed  and  a 
factor  of  2n;  they  correspond  to  those  given  by  Hocking  (1985).  This 
interpretation  assumes  the  Taylor  hypothesis,  i.e.,  the  eddies  causing  the 
diffraction  pattern  are  "frozen"  and  do  not  evolve  with  time  as  they  pass  over 
the  receiver  array..  Power  law  behavior  in  the  viscous  region  is  commonly 
described  by 
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o  E(k)<*  k‘7( 

(5.2) 

where  k  Is  the  horizontal  wave  number.  The  region  over  which  the  power  law 
behavior  is  valid  increases  with  height  because  so  does  the  length  scale  for 
the  viscous  region.  Note  that  the  slopes  for  the  power  law  behavior  (see  Table 
5.4)  for  76  and  79  km  are  almost  as  large  as  -7.  This  behavior  does  not  hold  for 
altitudes  above  79  km. 

The  region  between  approximately  10  and  200  s  is  the  inertial  region. 
This  region  exhibits  much  weaker  power  law  behavior  than  the  viscous 
region.  Some  power  law  curve  fits  for  this  region  are  shown  in  Figures  5.14  - 
5.17.  The  radar  does  directly  detect  eddies  on  the  length  scales  of  the  inertial 
region.  This,  combined  with  fluctuations  of  the  wind  speeds  on  time  scales 
equivalent  to  the  inertial  region  length  scales,  may  explain  the  departure  of 
the  power  law  curve  fit  from  the  commonly  expected  -5/3  behavior. 

( 

5.3  interpretation 

( 

The  first  step  in  interpreting  the  output  of  the  correlation  integral 
algorithm  lies  in  determining  the  slope  as  a  function  of  e.  Once  the  slope  has 
been  calculated,  we  can  determine  if  there  is  an  attractor  present  and  its 
f  dimension.  First,  let  us  discuss  how  to  interpret  the  slope  as  a  function  of  e. 

5.3.1  Interpreting  the  Slope  as  a  Function  of  e 

i  The  slope  of  the  lines  from  the  correlation  integral  algorithm  should 

equal  the  fractal  dimension  of  the  attractor  when  plotted  on  a  In-ln  graph.  The 
plot  of  the  slope  of  the  natural  logarithm  of  the  correlation  integral  against 

the  natural  logarithm  of  the  radius  e  can  be  divided  into  four  distinct  regions. 

f 


C 
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Each  region  yields  information  on  the  algorithm  and  the  underlying  attractor. 
Figure  5.18  shows  a  idealized  diagram  of  such  a  plot. 

Region  A  in  Figure  5.18  is  the  section  where  the  statistics  are  limited  by 
the  number  of  points  in  the  data  set.  The  slope  in  this  region  starts  out  near 
zero  and  often  is  characterized  by  some  oscillatory  behavior.  The  slope  is  zero 
because  the  number  of  points  in  a  "ball"  of  radius  e  approaches  zero  as  e 
approaches  zero.  Oscillatory  behavior  is  the  result  of  fluctuations  in  the 
counting  statistics  because  of  the  limited  number  of  data  points. 

Region  B  is  dominated  by  instrumental  and  external  noise.  The  slope  is 
proportional  to  the  embedding  dimension.  In  this  region,  the  "balls"  are 
smaller  than  the  smallest  temporal  and  spatial  scale  of  the  attractor.  Often,  the 
slope  will  only  approach  the  embedding  dimension;  this  occurs  primarily  at 
the  higher  embedding  dimensions  because  there  aren't  enough  points  to 
adequately  saturate  higher  dimensional  spaces. 

There  will  always  be  a  transition  between  regions  A  and  B  for  data  sets 
where  there  arc  a  finite  number  of  points.  As  the  number  of  points  increases 
this  transition  would  occur  at  increasingly  smaller  values  of  e. 

Region  C  is  the  "scaling"  region.  This  region  is  characterized  by 
constant  slope  which,  at  high  enough  embedding  dimension,  should  equal  the 
fractal  dimension  of  the  attractor.  The  slopes  should  converge  to  a  common 
value  for  embedding  dimensions  equal  or  greater  than  2d+l  since  this 
dimension  phase  space  will  completely  embed  the  attractor.  The  width  of  the 
scaling  region  is  vital.  The  limits  imposed  on  the  original  definition  of  the 
correlation  integral  necessitate  a  scaling  region  which  spans  several  orders  of 
magnitude.  Noise  reduces  the  width  of  the  scaling  region  and  may  make  it 
difficult  to  discern.  A  high  signal-to-noise  ratio  yields  a  larger  scaling  region 
whereas  signals  characterized  by  a  low  signal-to-noise  ratio  will  have  a  very 
small  scaling  region  or  may  not  have  one  at  all. 
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There  may  be  multiple  scaling  regions  present  in  region  C  if  there  are 
attractors  of  different  length  scales  present  in  the  data.  Multiple  attractors  in 
the  data  would  lead  to  distinct  plateaus  in  the  slope  diagram. 

Region  D  is  characterized  by  the  slope  converging  to  zero.  As  e 
increases  and  approaches  the  size  of  the  attractor,  the  "balls”  contain  nearly 
the  total  number  of  points  in  the  system.  Thus  the  slope  approaches  zero.  In 
some  systems,  the  transition  from  region  C  to  D  is  characterized  by  a  small 
"hump"  or  region  of  higher  slopes.  This  "hump"  occurs  because  of  edge  effects 
of  the  "balls"  used  to  count  the  points  and  curvature  of  the  attractor 
(Brandstater  and  Swinney,  1987). 

5.3.2  Calculating  the  Slope 

The  slope  of  the  lines  in  Figure  5.9  was  calculated  using  a  seven  point 
least  squares  fit  (the  algorithm  was  taken  from  Press  ct  al.,  1986).  The  slope 
calculated  for  each  point  was  the  least  squares  fit  of  the  point  extending  three 
points  on  either  side.  The  errors  of  the  least  squares  fit  were  calculated  but  not 
displayed  so  as  to  render  the  figures  more  legible. 

The  results  of  this  analysis  are  shown  in  Figure  5.19.  Examination  of  the 
figures  reveals  no  apparent  scaling  region  for  any  level.  The  data  set  seems  to 
be  characterized  by  noise  for  times  scales  less  than  6  min  40  s  (i.e.,  1500 
points).  There  may  be  a  higher  dimension  attractor  in  the  data  but  it  can  not 
be  detected  with  only  1500  points  and  thus  will  appear  as  noise.  While 
recognizing  this  possibility,  we  will  consider  the  results  to  be  noise  for 
purposes  of  discussion.  The  slopes  of  the  entire  correlation  integral  analysis 
are  given  in  Appendix  B. 
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5.3.3  Whv  Noise? 


There  does  not  appear  to  be  an  attractor  in  this  data  set.  Noise  seems  be 
present  in  the  signal  for  all  time  scales  that  were  investigated.  This  leads  us  to 
the  question  —  why  have  we  found  only  noise  and  should  we  have  expected  it? 

The  lack  of  positive  results  may  be  attributed  to  one  or  any  combination 
of  reasons.  These  fall  into  three  broad  categories:  difficulties  with  the 
algorithm,  difficulties  with  the  data  and  an  invalid  hypothesis. 

5.3.3. l  Difficulties . with  the  algorithm 

As  was  previously  described  in  Chapter  III,  the  Grassberger-Procaccia 
algorithm  has  a  number  of  weaknesses  that  make  it  difficult  to  implement. 
Among  these  are  the  finite  size  of  the  data  set,  noise  (both  instrumental  and 
external)  and  the  inability  of  the  algorithm  to  detect  the  presence  of  an 
attractor  amid  a  wide  spectrum  of  competing  signals. 

The  idealized  models  of  gravity  wave  breaking  invoke  a  single 
monochromatic  wave  that  becomes  unstable  at  some  point  as  it  travels  upward. 
This  is  almost  certainly  not  the  case  in  the  atmosphere;  instead  of  a  single 
monochromatic  gravity  wave,  there  is  an  entire  spectrum  of  gravity  waves 
which  travel  in  packets.  As  pointed  out  earlier,  the  spectrum  of  wr.ves  gives 
rise  to  nonlinear  wave-wave  interaction. 

It  is  not  known  whether  the  Grassberger-Procaccia  algorithm  could 
pick  out  the  presence  of  a  strange  attractor  in  such  a  sea  of  competing  signals, 
some  of  which  might  have  time  scales  on  the  order  of  the  attractor.  While 
there  has  been  a  limited  investigation  into  the  use  of  the  algorithm  in  the 
presence  of  two  competing  signals,  its  use  on  a  spectrum  of  signals  remains 
largely  unexplored. 

The  limited  size  of  the  data  set  could  also  prevent  detection  of  an 
attractor,  especially  if  the  dimension  is  greater  than  3.  The  1500  point  case  is 
not  adequate  to  accurately  estimate  the  dimension  of  systems  with  dimensions 
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greater  than  or  equal  to  2.  However,  1500  points  should  be  sufficient  to  test  for 
the  presence  of  an  attractor  of  dimension  3  or  less. 

We  did  perform  a  limited  set  of  calculations  using  5000  points  for  82  and 
85  km.  The  data  analyzed  was  an  average  of  all  four  antennas.  Unfortunately, 
these  calculations  were  performed  with  an  incorrect  value  of  the  time  delay 
used  to  reconstruct  the  phase  space  vectors;  the  value  of  the  time  delay  was 
only  2  instead  of  the  values  given  in  Table  5.1.  The  results  of  these  calculations 
are  shown  in  Figure  5.20.  From  the  plot  of  the  slope  as  a  function  of  e,  there 
doesn't  appear  to  be  an  attractor  in  this  data  set  for  the  5000  point  time  scale 
(i.e.,  22  min  13s).  While  it  is  true  that  these  results  are  not  as  reliable  as  those 
for  the  1500  point  cases  because  of  the  incorrect  choice  of  time  delay,  the 
greater  number  of  points  makes  the  exact  choice  of  the  time  delay  less  crucial, 
especially  since  we  found  only  noise.  Had  the  5000  point  case  suggested  the 
presence  of  an  attractor,  we  would  be  inclined  to  interpret  the  results  with 
more  caution.  We  look  upon  the  5000  point  case  as  further  support  for  the 
absence  of  an  attractor  and  there  only  being  noise  in  this  data  set. 

The  degree  of  autocorrelation  and  the  lack  of  general  independence  as 
determined  from  the  mutual  information  also  reduces  the  useable  size  of  the 
data  set  and  contributes  to  the  problems  in  the  algorithm  from  a  limited 
number  of  points. 

The  noise  in  the  signal  is  the  biggest  problem.  There  may  only  be  noise 
in  the  data  on  the  time  scales  examined  here  or  the  signal-to-noise  ratio  may 
be  so  low  as  to  completely  obscure  the  scaling  region  of  an  attractor.  The 
degree  to  which  the  noise  is  instrumental  or  external  will  be  touched  upon  in 
a  later  section. 

5. 3. 3.2  Difficulties. -  with  the  Data 

The  problems  associated  with  the  data  fall  into  two  categories:  no 
gravity  waves  in  the  data  set  and  the  quality  of  the  data  itself.  The  issue  of  data 
quality  overlaps  some  of  the  limitations  of  the  Grassberger-Procaccia 
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algorithm  that  have  been  discussed  in  previous  sections.  The  second  problem 
is  not  really  a  problem  but  a  statement  of  fact:  we  found  nothing  but  noise  in 
the  data  because  there  really  was  nothing  but  noise  in  the  data.  First  let  us 
discuss  the  quality  of  the  data  itself. 

The  biggest  problem  that  occurred  was  the  high  frequency  "ringing" 
present  in  the  data  in  the  second  and  third  file.  This  so  contaminated  the  data 
that  some  preliminary  analysis  of  these  signals  suggested  the  presence  of  an 
attractor  with  finite  dimensions  until  the  cause  was  discovered.  (The  plots  of 
the  original  data  set  that  were  supplied  to  us  indicated  no  such  rroblems).  As 
mentioned  earlier,  this  high  frequency  contamination  of  the  data  precluded 
the  use  of  almost  an  hour  of  data,  thus  severely  reducing  the  size  of  the  data 
set.  This  exacerbated  the  problems  that  occur  with  finite  data  sets  as  described 
in  previous  sections. 

addition  to  the  "ringing",  gain  changes  and  drift  in  the  receiver 
presented  an  additional  source  of  error.  While  each  of  these  problems  could 
have  been  corrected,  doing  so  would  have  introduced  an  additional  degree  of 
uncertainty  into  any  dimension  estimates.  These  problems  precluded  use  of  the 
last  two  files. 

The  elimination  of  almost  two  thirds  of  the  data  set  emphasizes  the 
Grassberger-Procaccia  algorithm's  need  for  lengthy  data  records.  Data  analysts 
always  complain  they  could  use  more  data,  but  in  this  case  it  seems  to  be  an 
absolute  necessity.  Without  a  sufficient  number  of  points,  the  entire  algorithm 
is  not  statistically  valid. 

This  data  set  may  contain  no  evidence  of  gravity  wave  breaking.  The 
data  set  was  originally  suggested  to  us  as  one  that  potentially  contained 
evidence  of  gravity  wave  breaking. 

Meek  and  Manson  (1987)  used  this  data  set  to  examine  "glints"  in  the 
pattern  of  echoes  at  the  receiver.  The  "glints"  were  regions  of  strong 
persistent  echoes  that  travelled  in  a  straight  line  through  the  receiver  array. 


They  speculate  that  the  glints  were  either  long  lived  pockets  of  turbulence  or 
focusing  of  the  signal  due  to  favorable  alignment  of  surfaces  of  the  wave. 

Meek  and  Manson  (1987)  found  two  glints  in  the  data  set  from  day  214 
that  has  been  examined  in  this  study.  The  first  glint  was  observed  between 
18:31:24  and  18:34:36  GMT  as  it  travelled  through  the  receiver  array.  The 
second  glint  was  observed  between  18:47:48  and,  18:49:24  GMT.  The  velocity  of  a 
glint  could  be  deduced  by  tracking  the  amount  of  time  it  took  to  cross  the 
receiver  array;  Meek  and  Manson  (1987)  estimated  that  both  glints  had 
velocities  of  ~  50  m  s'1.  If  the  glints  were  caused  by  focusing  of  the  turbulence 
along  the  surface  of  a  wave,  the  period  of  the  wave  can  be  deduced  from  the 

difference  in  time  between  the  glints.  The  period  of  a  wave  travelling  50  m  s*1 

is  approximately  5  minutes  which  is  roughly  the  same  as  estimates  of  the 
Brunt-VaisSUa  period  in  the  mesosphere. 

Meek  and  Manson  (1987)  used  a  wave  model  to  determine  if  the  glints 
could  be  caused  by  focusing  along  the  surface  of  a  travelling  wave.  The  wave 
model  reproduced  the  observations  of  the  glints  well.  Meek  and  Manson 

further  speculated  that  the  waves  responsible  for  the  glints  were  close  to  their 

critical  levels  since  their  phase  velocities  matched  the  parallel  component  of 
the  background  wind  (an  instability  criterion  given  by  Fritts  and  Rastogi, 
1985).  Thus,  while  Meek  and  Manson  did  not  directly  observe  gravity  waves 
breaking,  they  speculated  the  conditions  were  favorable  for  this  occurrence. 

Unfortunately,  we  could  find  no  evidence  of  gravity  wave  activity  in 
our  analysis  of  the  data  set.  Power  spectra  of  the  data  for  each  level  and 
antenna  showed  no  evidence  for  wave  activity  near  the  Brunt-VSisSlia 
frequency.  To  detect  periods  in  the  power  spectrum  out  to  the  Brunt-VSisalia 
period  required  using  data  from  the  contaminated  second  and  third  files.  The 
problems  previously  discussed  in  detail  may  h3vc  prevented  us  from  detecting 
long  period  signals  in  the  power  spectrum. 

Meek  and  Manson  (1987)  found  evidence  of  wave  activity  in  this  data 
set.  Furthermore,  they  speculated  that  the  waves  causing  the  glints  may  be 
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near  saturation.  While  our  analysis  of  the  data  set  could  not  support  this 
conclusion  for  the  reasons  described  above,  neither  could  we  rule  it  out 
conclusively.:  This  still  leaves  us  with  the  unanswered  question  of  whether  or 
not  the  presence  of  breaking  gravity  waves  exists  in  the  data  set. 

5.3.3. 3  Invalid  Hypothesis 

The  last  reason  why  we  may  not  have  detected  anything  but  noise  in  the 
signal  is  that  the  hypothesis  is  invalid.  Gravity  wave  breaking  and  the 
subsequent  decay  to  turbulence  may  not  be  characterized  by  chaotic  behavior. 
The  analogies  to  laboratory  studies  may  be  incorrect  or  they  may  not  extend  to 
the  real  atmosphere  where  there  is  an  entire  spectrum  of  activity  besides 
gravity  waves.  On  the  other  hand,  the  hypothesis  may  be  valid;  the  breakdown 
of  atmospheric  gravity  waves  may  be  characterized  by  chaotic  behavior  but  it 
is  of  such  an  ephemeral  nature  we  may  not  be  able  to  detect  it. 

The  extension  of  the  results  from  laboratory  experiments  to  the 
atmosphere  may  be  invalid.  In  our  review  of  some  of  the  laboratory  studies,  we 
examined  chaotic  behavior  which  arose  from  convective  instabilities 
(Rayleigh-Benard  convection)  and  dynamic  instabilities  (Couette-Taylor  flow) 
in  closed  systems  and  dynamic  instabilities  in  open  systems  (the  excited  jet). 
None  of  these  regimes  are  similar  to  the-  type  of  flow  and  instabilities  thought 
to  describe  gravity  waves  and  their  breakdown,  although  the  shear 
instabilities  present  in  the  excited  jet  are  somewhat  similar  to  one  type  of 
dynamical  instability  that  may  occur  in  gravity  waves.  It  may  be  wrong  to 
infer  chaotic  behavior  in  breaking  gravity  waves  on  this  basis. 

Laboratory  experiments  are  held  under  tightly  controlled  conditions 
which  are  very  unlike  those  found  in  the  atmosphere.  Only  one  spectrum  of 
activity  is  studied  making  it  easier  to  detect  chaotic  behavior  in  a  fluid.  Gravity 
wave  breaking  in  the  atmosphere  is  accompanied  by  a  wide  range  of  other 
types  of  processes  in  the  atmosphere,  some  of  which  have  time  scales  that 
overlap  those  of  gravity  waves.  Detecting  chaotic  behavior  in  such  a  weher  of 


potentially  conflicting  data  may  be  nearly  impossible.  Furthermore,  the 
interaction  of  gravity  wave  breaking  and  other  types  of  atmospheric  motion 
may  eliminate  any  chaotic  behavior  that  might  occur  otherwise. 

Laboratory  experiments  allow  for  the  generation  of  a  long  period  of 
measurements.  The  conditions  which  control  the  instability  of  the  fluid  flow 
can  be  held  constant  allowing  measurements  to  be  made  of  the  long  term 
evolution  of  its  chaotic  behavior.  The  flow  is  "continuously  unstable". 

Gravity  wave  breaking  is  a  transient  phenomenon;  gravity  waves 
propagate  upward  and  may  become  unstable  and  break  down.  \  better 
extension  of  the  laboratory  analogy  would  be  a  continuous  source  of  gravity 
waves  which  become  unstable  and  break  without  modifying  the  basic  state 
flow  at  some  level.  This  scenario  is  unlikely  to  ever  occur  in  the  real 
atmosphere,  much  less  at  a  time  and  place  where  measurements  were  being 
made. 

Chaotic  behavior  in  gravity  wave  breaking  may  be  so  ephemeral  that 
we  may  never  be  able  to  detect  it.  Bonetti  and  Boon  (1989)  recommended  a 
sampling  rate  of  10  to  30  measurements  per  pseudo-period  of  the  orbit  of  the 
attractor.  If  we  consider  Smith's  (1988)  estimate  of  the  number  of  points 
required  to  accurately  estimate  the  dimension  using  the  Grassberger- 
Procaccia  algorithm  for  an  attractor  of  dimension  2,  and  Bonetti  and  Boon’s 
estimate  of  the  sampling  rate,  we  must  sample  approximately  between  59  and 
176  orbits  to  fully  characterize  the  attractor!  If  we  assume  that  the  pseudo¬ 
period  of  an  attractor  in  the  breaking  of  gravity  waves  is  on  the  order  of  the 
Brunt-Vaisiilla  period  (5  minutes)  and  a  dimension  on  the  order  of  2,  we  would 
have  to  make  measurements  for  more  than  6  hours..  The  length  of  time  may  be 
even  greater  if  the  attractor  were  of  higher  dimensions  which  is  very  likely. 

Six  hours  is  practically  an  eternity  for  the  phenomena  we  are 
considering  in  the  atmosphere.  The  data  set  would  not  be  stationary  over  this 
length  of  time  and  atmospheric  flow  regimes  with  longer  time  scales  (tides, 
synoptic  scale  activity)  would  contaminate  the  data  set.  This  requirement 
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almost  certainly  eliminates  finding  chaotic  behavior  in  gravity  wave 
breaking  since  gravity  waves  that  break  decay  to  turbulence  at  a  rate  much 
faster  than  6  hours. 
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Figure  5.1  Mutual  information  for  antenna  #1  calculated  from  the  first 
6000  points  for  levels  76  through  94  km. 
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Figure  5.2  Mutual  information  for  antenna  #2  calculated  from  the  first 
6000  points  for  levels  76  through  94  km.. 


Figure  5.4  Mutual  information  for  antenna  #4  calculated  from  the  first 
6000  points  for  levels  76  through  94  km. 


Figure  5.5  (c) 
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Autocorrelation 


Figure  5.6(a) 


Figure  5.6  Autocorrelation  (2000  point  groups)  for  antenna  #2  for  (a)  76 
km,  (b)  79  km,  (c)  82  km,  (d)  85  km,  (e)  88  km,  (0  91  km  and  (g) 
94  km. 
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Figure  5.7  (d) 


Lag 


Figure  5.7  (e) 
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Figure  5.8(a) 


Figure  5.8  Autocorrelation  (2000  point  groups)  for  antenna  #4  for  (a)  76 
km,  (b)  79  km,  (c)  82  km,  (d)  85  km,  (e)  88  km,  (0  91  km  and  (g) 
94  km. 
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Figure  5.8  (g) 
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Figure  5. 


Figure  5.10(a) 


Power  spectral  density  for  antenna  #1  for  (a)  76  km,  (b)  79  km, 


Figure  5.10  (d) 
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Figure  5.10  (g) 
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Figure  5.12(b) 
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Figure  5.12  (c) 
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Figure  5.13  (c) 
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Power  Spectral  Density  Power  Spectral  Density 
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Figure  5.13  (f) 
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Figure  5.13  (g) 
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Figure  5.14  (c) 
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Figure  5.14  (h) 
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Figure  5.15  (a) 
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Figure  5.15(b) 

Figure  5.15  Power  law  fits  for  the  power  spectra  from  antenna  #2  for  (a)  76 
km,  (b)  76  km,  (c)  79  km,  (d)  79  km,  (e)  82  km,  (0  82  km,  (g)  85 
km,  (h)  85  km,  (i)  88  km,  (j)  91  km  and  (k)  94  km.. 
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Figure  5.17  (i) 
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Figure  5.18  Idealized  depiction  of  the  slope  of  the  natural  logarithm  of  the 
correlation  integral  as  a  function  of  e.  The  figure  is  divided  into 
four  distinct  regions  as  indicated  by  the  dashed  lines. 


Figure  5.19  Slope  of  the  correlation  integral  plotted  as  a  function  of 

embedding  dimension  for  points  1  through  1500  from  antenna  1 
for  (a)  76  km;  (b)  79  km;  (c)  82  km;  (d)  85  km  (e)  88km;  (f)  91  km 
and  (g)  94  km. 
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Figure  5.19(b) 
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CHAPTER  VI 


CONCLUSION  AND  RECOMMENDATIONS  FOR  FUTURE  WORK 

6.1  Conclusions 

Wc  did  not  detect  a  strange  attractor  with  dimension  less  than  three  in 
the  data  from  the  Saskatoon  partial  reflection  radar  for  the  time  scales  which 
were  studied.  While  we  can  not  assert  that  there  was  only  noise  in  the  data  set 
because  of  the  small  number  of  points  which  were  examined,  the  supporting 
evidence  from  the  power  spectra  suggest  that  we  investigated  time  scales 
primarily  in  the  viscous  and  inertial  region  and  that  the  dimension  of  this 
system  was  greater  than  3. 

The  limited  size  of  the  data  set  was  one  of  the  major  reasons  we  can 
definitively  rule  out  the  presence  of  an  attractor  of  only  less  than  dimension 
3.  All  indications  are  that  the  dominant  signal  was  due  to  noise.  We  interpret 
this  noise  to  be  isotropic  three  dimensional  turbuljnce. 

The  difficulties  in  detecting  a  strange  attractor  in  the  saturation  of 
middle  atmosphere  gravity  waves  have  been  discussed  in  Chapter  V.  It  does 
seem  likely  that  the  breaking  of  a  gravity  wave  (or  a  spectrum  of  gravity 
waves)  may  be  so  ephemeral  that  limitations  in  the  current  implementation  of 
the  Grassberger-Procaccia  algorithm  will  prevent  us  from  ever  detecting  it. 
The  decay  to  turbulence  should  occur  quickly  enough,  even  in  the  slower 
slantwise  static  instability  mechanism,  that  there  will  not  be  a  statistically 
significant  number  of  orbits  of  the  attractor  to  be  useful.  At  least  50  orbits  of 
the  attractor  are  necessary  to  implement  the  Grassberger-Procaccia  algorithm. 
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An  even  bigger  consideration  remains  -  if  the  attractor  is  so  ephemeral  is  it 
worth  the  trouble  to  study? 

A  similar  question  remains  as  to  what  a  strange  attractor  might  mean 
physically.  In  the  case  of  the  absorption  of  gravity  waves,  it  still  may  provide 
a  clue  as  to  a  mechanism  by  which  the  wave  becomes  unstable.  As  pointed  out 
in  previous  discussions,  the  currently  accepted  conceptual  models  for  gravity 
wave  breaking  are  most  likely  wrong  if  not  merely  inaccurate.  Detecting  the 
presence  of  a  strange  attractor  and  its  dimension  would  provide  insight  into 
the  route  to  chaos  in  this  system  and  perhaps  insight  into  the  nature  of  the 
mechanism  of  gravity  wave  saturation. 

If  there  is  an  attractor  associated  with  the  nonlinear  wave-wave 
interaction,  does  this  reveal  anything  about  the  mechanism  behind  gravity 
wave  saturation  or  just  something  about  the  particular  spectrum  of  waves 
undergoing  saturation  in  the  data  set?  We  must  be  cautious  in  generalizing 
any  future  results  to  gravity  wave  saturation  itself. 

6.2  Recommendations  For  Future....  Work 

This  study  raises  many  more  questions  than  it  has  definitively 
answered.  While  it  seems  unlikely  that  an  attractor  will  be  found  which 
describes  the  transition  to  chaos  over  a  short  time  scale  in  meteorological  data, 
more  work  is  needed  before  it  can  be  absolutely  discounted.  The  suggestions 
for  further  study  fall  into  three  areas:  work  on  the  technique,  extension  of 
this  work  with  better  radar  data  and  extension  of  this  work  with  boundary 
layer  data. 

The  Grassberger-Procaccia  algorithm  and  its  extensions  have  been  used 
successfully  in  a  number  of  different  theoretical  and  laboratory  studies  of 
chaotic  behavior.  Its  weaknesses  have  been  well  studied  and  documented  in 
the  literature,  for  the  most  part.  We  are  uncertain  whether  this  algorithm 
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could  detect  the  presence  of  a  strange  attractor  in  the  confusing  welter  of 
competing  signals  that  make  up  a  typical  meteorological  data  set.  The  degree  of 
sensitivity  of  the  Grassberger-Procaccia  algorithm  to  a  signal  that 
encompasses  a  wide  spectrum  of  different  types  of  behavior  (periodic,  noise 
and  chaotic)  must  be  established. 

The  Grassberger-Procaccia  algorithm  must  be  investigated  as  to  its 
ability  to  detect  a  strange  attractor  amidst  noise  and  a  spectrum  of  periodic 
signals  such  as  are  found  in  a  meteorological  data  set.  This  investigation  must 
determine  the  relative  strength  of  the  chaotic  signal  that  is  necessary  to  be 
detected.  It  also  must  determine  the  number  of  points  needed  to  make  it  work. 

This  study  should  be  extended  to  other  middle  atmosphere  data  sets  so 
that  its  conclusions  may  be  confirmed.  Such  data  sets  must  consist  of  much 
longer  records  without  the  gain  changes  and  contamination  of  the  signal  that 
was  found  in  the  data  used  in  this  study.  It  would  also  be  useful  to  employ  this 
technique  on  the  radar  derived  winds  which  are  closer  to  a  truer  description 
of  gravity  waves.  The  nebulous  connection  between  the  signals  that  are 
measured  by  partial  reflection  radars  and  the  middle  atmosphere  may  obscure 
the  presence  of  chaotic  behavior.  Wind  measurements  deduced  from  the  radar 
signal  do  carry  a  degree  of  smoothing  that  could  complicate  the  analysis,  but 
in  general  they  will  be  more  representative  uf  gravity  wave  behavior. 

Other  algorithms  should  be  used  to  supplement  the  Grassberger- 
Procaccia  algorithm.  The  improved  box  counting  algorithm  (Liebovitch  and 
Toth,  1989)  should  be  examined  to  determine  if  it  is  a  useful  adjunct  to  the 
Grassberger-Procaccia  algorithm.  The  nearest  neighbor  method  of  Badii  and 
Politi  (1987)  should  be  tried  as  well  as  some  of  the  additional  extensions  of  the 
Grassberger-Procaccia  algorithm  (Franaszek,  1989;  Ellner,  1988).  Higuchi 
(1988)  has  suggested  an  algorithm  using  a  fractal  length  of  curve  technique 
which  shows  promise  for  low  (<2)  dimension  attractors;  this  might  be  tried  as 
well. 
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The  atmospheric  boundary  layer  would  be  a  logical  place  in  which  to 
extend  this  study.  Like  the  middle  atmosphere,  it  is  the  home  to  unstable  wave 
activity  that  decays  to  turbulence  (e.g.,  Kelvin-Helmholtz  waves).  The  study  of 
turbulence  in  the  boundary  layer  is  also  well  advanced.  Most  importantly, 
measurements  of  the  boundary  layer  are  of  a  much  higher  quality  than  those 
of  the  middle  atmosphere.  Most  boundary  layer  parameters  can  be  sampled  in 
situ  rather  than  deduced  by  remote  sensing  techniques  as  is  commonly  done  in 
the  middle  atmosphere.  The  amount  of  external  noise  in  boundary  layer  data  is 
also  much  less  than  that  for  the  middle  atmosphere.  Measurements  of  several 
different  variables  can  be  made  which  help  provide  a  better  understanding  of 
the  dynamics  of  an  attractor  if  one  is  found.  In  addition,  measurements  of  the 
boundary  layer  can  be  made  at  much  higher  sampling  rates  and  for  longer 
continuous  periods  than  in  the  middle  atmosphere. 

Stationary  data  sets  are  perhaps  an  even  greater  problem  in  the 
boundary  layer  than  in  the  middle  atmosphere.  The  boundary  layer  undergoes 
a  tremendous  change  during  the  diurnal  heating  cycle,  exhibiting 
remarkably  different  types  of  behavior  between  daytime  and  nighttime.  The 
boundary  layer  study  of  Tsonis  and  Eisner  (1988)  was  flawed  not  only  because 
it  did  not  consider  enough  points  in  applying  the  Grassberger-Procaccia 
algorithm  and  the  data  set  was  highly  non-stationary,  but  because  it  offered  no 
reason  (i.e.,  physical  insight)  as  to  why  there  might  be  a  strange  attractor  in 
meteorological  data  over  a  short  time  scale. 

There  may  a  greater  chance  of  finding  chaotic  behavior  in  the 
structure  of  turbulence  rather  than  the  transition  from  laminar  to  turbulent 
flow.  Recent  work  using  multi-fractals  suggests  that  the  structure  of 
turbulence  is  fractal  in  nature  (Chhabra  et  al.„  1989;  Meneveau  and  Nelkin, 
1989;  Smith  et  al..  1986).  There  hints  of  this  in  the  work  here  for  time  scales  in 
the  viscous  region.  This  is  another  avenue  for  research  that  should  be 
explored. 
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APPENDIX  A 


Appendix  A  contains  the  figures  depicting  the  results  of  the  correlation 
integral  calculations  described  in  Chapter  V  for  all  antennas  and  levels.  Output 
from  only  every  other  embedding  dimension  are  shown  in  an  effort  to  render 
the  figures  more  legible. 
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APPENDIX  B 


This  appendix  contains  figures  depicting  the  slopes  of  the  correlation 
integral  calculations  shown  in  Appendix  A.  The  slopes  were  calculated  using  a 
seven  point  fit  as  described  in  Chapter  V.  Output  from  only  every  other 
embedding  dimension  are  shown  in  an  effort  to  render  the  figures  more 
legible. 
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